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In the filter evaluation process, post - filtered pulse-pair 
width estimates and power levels are also used to measure the 
effectiveness of the filters. 

The results presented support the use of an adaptive 
clutter rejection filter for reducing the clutter induced 
bias in pulse-pair estimates of windspeed. The adaptive 
clutter rejection filter is also shown to perform better than 
one clutter rejection filter commonly used in land based air 
traffic control radar systems. 



TABLE OF CONTENTS 


Page 


TITLE PAGE i 

ABSTRACT ii 

ACKNOWLEDGEMENTS iv 

LIST OF TABLES viii 

LIST OF FIGURES ix 

CHAPTER 


I. MICROBURST DETECTION 1 

Microburst 1 

Pulsed Doppler Weather Radar 4 

Pulse-Pair Estimator 9 

F-Factor 12 

Statement of the Problem 13 

II. CLUTTER MODELING 16 

Introduction 16 

ARMA Modeling 18 

AR Modeling 19 

Levinson-Durbin Algorithm 20 

Linear Prediction Error Filter 21 

Lattice Prediction Error Filter 24 

Coefficient Determination 27 

Burg ' s Block Implementation 27 

Adaptive Implementation 28 

Gradient Adaptive 28 

Least Squares Adaptive 3 3 

Model Order 38 

III. THE EFFECTS OF CLUTTER IN ESTIMATING 

WINDSPEED 39 

NASA Model 39 

Bias in Pulse -Pair Mean Estimates Due to 

Clutter 40 

Analysis of a Dry Microburst 42 

Validation of Mean Estimates 4 5 

Width Estimates 4 5 

Power Levels 4 5 


VI 


Table of Contents (Continued) 

Page 

Threshold Levels 48 

Clutter Analysis 49 

IV. CLUTTER REJECTION FILTERS 54 

The Optimum Clutter Rejection Filter 54 

The Complex Square Root Normalized 
Recursive Least Squares Lattice 

Estimation Algorithm 55 

10th Order AR Clutter Model 60 

10th Order Complex Coefficient 

FIR Filters 63 

Analysis of Filtering Schemes 65 

Filters Based on the Clutter 

in Each Range Cell 65 

Filter Based on the Clutter 

in a Single Range Cell 69 

Comparison with a Pulse Canceller 75 

V. CONCLUSIONS 81 

APPENDICES 84 

A. The Complex Square Root Normalized 

Recursive Least Squares Lattice 

Estimation Algorithm 85 

B. NASA Model Parameters 127 

C. Fourier Spectral Estimates of a Dry 

Microburst Plus Clutter 129 

D. Fourier Spectral Estimates of a Dry 

Microburst Without Clutter 135 

E. Fourier Spectral Estimates of Clutter 141 

F. Magnitude Response of the 10th Order 

FIR Filters 147 

G. Phase Response of the 10th Order 

FIR Filters 153 

H. Filtered Spectrum Using the Appropriate 

Model Based Filter Coefficients 159 



VI 1 


Table of Contents (Continued) 

Page 

I. Filtered Spectrum Using the Filter 

Coefficients for Range Cell 20 165 

J. Filtered Spectrum Using a Pulse 

Canceller Filter 171 


LITERATURE CITED 


177 



LIST OF TABLES 

Table Page 

I. Levinson-Durbin Algorithm 22 

II. Gradient Adaptive Lattice Algorithm 30 

III. Algorithm for Determining the AR Model 
Coefficients Based on a Gradient 

Adaptive Lattice Structure 31 

IV. Complex Least Squares Adaptive 

Lattice Algorithm . . , , 34 

V. Algorithm for Determining the AR Model 
Coefficients Based on a Complex LS 

Adaptive Lattice Structure 36 

VI. Complex Square Root Normalized Recursive 

Least Squares Lattice Estimation 

Algorithm 56 

VII. Algorithm for Determining the AR Model 

Coefficients Based on a Normalized 
Complex LS Adaptive Lattice 

Structure 61 

A- 1. Complex Least Squares Adaptive 

Lattice Algorithm 104 

A- II. Complex Square Root Normalized Recursive 

Least Squares Lattice Estimation 

Algorithm 113 

A- III. Algorithm for Determining the AR Model 

Coefficients Based on a Complex LS 

Adaptive Lattice Structure 121 

A- IV. Algorithm for Determining the AR Model 

Coefficients Based on a Normalized 
Complex LS Adaptive Lattice 
Structure 


125 



LIST OF FIGURES 


Figure Page 

1. Illustration of the outflow from a 

microburst 2 

2. "S" curve associated with a microburst 5 

3. Doppler return from a moving target 6 

4. Pulse modulation for a Doppler radar 8 

5. Block diagram of a coherent pulsed 

Doppler radar 10 

6. Fourier spectral estimate of a range cell 

containing microburst and clutter 

return . 17 

7. Linear prediction error filter implemented as 

a tapped delay line 23 

8. Linear prediction error filter implemented as 

a lattice structure 26 

9 . Transfer function implemented as a lattice 

structure for the gradient adaptive 

prediction error filter case 32 

10. Least squares adaptive lattice prediction 

error filter 35 

11. Transfer function implemented as a lattice 

structure for the least squares adaptive 

prediction error filter case 37 

12. Pulse-pair mean estimates of the return from 

a wet microburst and clutter 41 

13. Pulse-pair mean estimates of the return from 

a dry microburst and clutter 41 

14. Pulse-pair mean estimates of the return from 

a dry microburst 43 

15. Fourier spectral analysis of range cells 25, 

27, and 29 containing a dry microburst 44 


X 


List of Figures (Continued) 

Page 

16. Fourier spectral estimate of range cell 15 

with a low SNR 4 6 

17. Pulse-pair width estimates of the return from 

a dry microburst 47 

18. Power levels associated with the return from 

a dry microburst 47 

19. Pulse-pair mean estimates of the return from 

clutter 50 

20. Fourier spectral estimates of range cells 6, 

20, and 35 containing clutter 51 

21. Pulse-pair width estimates of the return from 

clutter 53 

22. Power levels associated with the return from 

clutter 53 

23. Complex square root normalized least squares 

adaptive lattice prediction error filter 57 

24. Transfer function implemented as a lattice 

structure for the complex square root 
normalized least squares adaptive 

prediction error filter case 59 

25. Fourier spectral estimate of the clutter in 

range cell 6 62 

26. AR spectral estimate of the clutter in range 

cell 6 62 

27. FIR clutter rejection filter implemented as 

a tapped delay line 64 

28. Pulse-pair mean estimates of filtered dry 

microburst plus clutter data where the 
appropriate filter based on the 
clutter in each range cell 
has been applied 


66 



XI 


List of Figures (Continued) 

Page 

29. Pulse-pair mean estimates of filtered dry 

microburst plus clutter data where the 

appropriate filter based on the 

clutter in each range cell has 

been applied compared to 

estimates given in the 

case of a dry 

microburst 

without 

clutter 66 

30. Filtered spectrum of range cell 7 68 

31. Pulse-pair width estimates after filtering using 

the filter coefficients obtained from 

the clutter in each range cell 68 

32. Power levels after filtering using the filter 

coefficients obtained from the clutter 

within each range cell 7 0 

33. Magnitude response of the filter designed for 

range cell 20 72 

34. Phase response of the filter designed for 

range cell 20 72 

35. Pulse-pair mean estimates after filtering 

using the filter coefficients for 

range cell 20 73 

36. Pulse-pair mean estimates after filtering using 

the filter coefficients for range cell 20 
compared to estimates given in the 
case of a dry microburst without 

clutter 73 

37. Pulse-pair width estimates after filtering using 

the filter coefficients for range 

cell 20 74 

38. Power levels after filtering using the filter 

coefficients for range cell 20 74 

39. Magnitude response of a pulse canceller 

filter 76 

40. Phase response of a pulse canceller filter 76 


Xll 


List of Figures (Continued) 

Page 

41. Pulse-pair mean estimates after filtering using 

a pulse canceller 77 

42. Pulse-pair mean estimates after filtering using 

a pulse canceller compared to estimates given 

in the case of a dry microburst without 

clutter 77 

43. Spectrum of range cell 24 after filtering 78 

44. Power levels after filtering using a pulse 

canceller 80 

45. Pulse-pair width estimates after filtering 

using a pulse canceller 80 

A-l. Decomposition of projections by oblique 

projection 95 

A-2. Least squares adaptive lattice prediction 

error filter 105 

A- 3. Complex square root normalized least squares 

adaptive lattice prediction error 

filter 114 

A- 4. Transfer function implemented as a lattice 

structure for the least squares adaptive 
prediction error filter case 122 

A- 5. Transfer function implemented as a lattice 

structure for the complex square root 
normalized least squares adaptive 
prediction error filter case 


126 



CHAPTER I 


MICROBURST DETECTION 
Microburst 

The contribution of microbursts to aircraft accidents 
was first recognized by Fujita [1] in the late 1970's while 
investigating the crash of Eastern 66 on June 24, 1975 at 
John F. Kennedy Airport in New York City. Between 1964 and 
1985, an estimated 26 major aircraft accidents resulting in 
626 fatalities and over 200 injures have been attributed to 
the microburst described by Fujita. Based on these figures, 
the FAA and NASA have proposed the use of airborne Doppler 
weather radar for the detection of a microburst in the 
vicinity of airport runways [2] . 

The term "microburst" was first used by Fujita to 
describe the vertical flow of wind which results in an 
horizontal outflow in all directions upon impact with the 
ground (Figure 1) . The horizontal outflow may extend in all 
directions within a 4km radius. The duration of a microburst 
is usually less than 10 minutes. This limited duration and 
the geographical boundaries of the horizontal outflow permit 
microbursts to go undetected by non-Doppler radars or 
ground-based anemometers [1] . Horizontal outflows extending 
in radius greater than 4km are termed macroburst and are 
easier to detect due to their large physical size. 




Figure 1. Illustration of the outflow from a microburst. 
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The aviational hazard associated with a microburst 
exists on final approach to a runway or on takeoff (Figure 
1). Headwinds are encountered by an aircraft on final 
approach as it moves through one side of the microburst. The 
pilot must then compensate to bring the aircraft down to the 
glide slope. At this point, the pilot has reduced his speed 
and is vulnerable to the unexpected downflow and tailwinds 
associated with the other side of the microburst. On 
takeoff, the downflow and tailwinds associated with the 
microburst combine to reduce the necessary lift required by 
the aircraft. 

Microburst can be subdivided into two categories based 
on the amount of precipitation present. A microburst 
containing high levels of precipitation is defined as being a 
"wet microburst". A low level of precipitation within the 
microburst leads to the term "dry microburst" . A high cloud 
base may allow time for evaporation in dry regions leading to 
the low levels present in the dry microburst. 

The precipitation present in a microburst provides the 
necessary targets for a Doppler weather radar which can be 
used to measure horizontal components of wind speed with 
respect to the aircraft. A wet microburst contains a 
sufficient target base to provide radar returns 
representative of the relative windspeed within the 
microburst. The detection of windspeed within a dry 
microburst is hampered due to the decreased target 


concentration . 
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A signature based on average horizontal components of 
windspeed relative to position within a microburst can be 
defined. This signature develops as the horizontal component 
of windspeed increases with distance from the center of the 
microburst. Figure 2 gives an example of the characteristic 
"S" curve which relates the average horizontal component of 
windspeed within a microburst to distance from the center of 
the microburst. The Doppler weather radar proposed by the 
FAA and NASA is intended to obtain information about average 
windspeed versus range and thus identify the presence of a 
microburst . 


Pulsed Doppler Weather Radar 

A Doppler weather radar is based on the Doppler effect 
defined by Christian Johann Doppler. The Doppler effect, 
applied to radar, relates the frequency shift in an 
electromagnetic signal to the relative motion of the target. 
The linear approximation between the frequency shift, f d , and 
the speed of the target, v, is given by 

f d = 2f t v/c. (1) 

The slope of this line is twice the ratio of the transmitted 
frequency, f j- , to the speed of light, c. Figure 3 
illustrates the shift in frequency of an electromagnetic 
signal incident upon a target with velocity, v. 

A Doppler radar can also be used to supply ranging 
information. This is accomplished by pulse modulating the 
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transmitted signal 


ft - f d 



negative Doppler shift 


positive Doppler shift 




Figure 3. Doppler return from a moving target. 
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transmitted waveform and gating the return signal at equal 
time intervals corresponding to successive target ranges 
(range bins or range cells) before the next pulse is 
transmitted. This is illustrated in Figure 4. The pulse 
repetition frequency, PRF, limits the maximum unambiguous 
range, R u> where 

R u = C/2PRF (2) 

shows that the transmitted signal must traverse its path to 
and from the target during the time interval allotted by the 
PRF. The ability to resolve two targets, a distance, D, 
apart, is also related to the PRF through 

D = c (t/ 2) (3) 
where r is the pulse width. 

The pulse modulation of the transmitted signal requires 
a coherent phase between pulses in order to extract frequency 
shift information. This is accomplished by maintaining a 
coherent oscillator (COHO) within the radar. The output from 
the COHO is at an intermediate frequency, IF, and is mixed 
with a stable local oscillator (STALO) to reach the required 
transmission frequency (RF) . The signal to be transmitted is 
then pulse modulated and amplified. 

The return signal, received during the gating process, 
is demodulated using the STALO in order to convert to the IF. 
In order to resolve the sign of the frequency shift, a second 
IF signal is formed which has a 90 degree phase shift with 
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(b) Gated return. 

Figure 4. Pulse modulation for a Doppler radar. 
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respect to the received IF signal. The two channels are 
termed inphase (I channel) and quadrature (Q channel) 
corresponding to the received IF signal and the 90 degree 
phase shifted version, respectively. A second demodulation 
is required to convert the IF signal to the Doppler frequency 
range. This is accomplished by mixing the COHO with the I 
and Q signals. The gated return is finally sampled and its 
value stored for future processing. Figure 5 shows a block 
diagram of a coherent pulsed Doppler radar [3] . 

The samples taken for a particular range cell over 
successive pulses are stored and then processed to determine 
the average target speed at that range. The pulse repetition 
frequency is also the sampling frequency for the range cell 
data since one sample is obtained per pulse per range cell. 
This sampling frequency, PRF, is small compared to the 
sampling frequency between successive range cells and 
therefore provides the smaller bandwidth needed to resolve 
the Doppler frequencies using conventional Fourier analysis 
or parametric spectral estimation. 

Pulse- Pair Estimator 

The precipitation present in a microburst yields 
multiple point targets within a range cell. The radar return 
from a range cell therefore consists of the sum of the 
returns from the individual targets (raindroplets) . Fourier 
spectral analysis of this range cell data reveals that the 




Figure 5. Block diagram of a coherent pulsed Doppler radar 
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estimated windspeed distribution can frequently be described 
as Gaussian in form with an associated mean and width [4] . 

The pulse -pair algorithm developed by Rummler [5] can be 
used to estimate the windspeed mean and width within a range 
cell. The algorithm uses estimates of the zeroth and first 
lag of the autocorrelation to determine the mean and width 
estimates. The pulse-pair width estimate, w, for the range 
cell data described previously can be written as 


w 


(PRF) 2 


2 

2 (n) 



I R (1) I 
R (0) 


(4) 


The autocorrelation lags in Equation 4 are estimated using 


1 N-l 

R (1) = - Z V* (n) V (n+1) (5) 

N n=0 


and 


1 N-l 

R (0 ) = - Z V* (n)V(n) (6) 

N n=0 

where V(n) is the sampled complex time series and N is the 
total number of samples. The pulse -pair mean estimate, f, 
takes the form 


(PRF) {ARG [R ( 1) ] } 
2n 


In comparing the pulse-pair method to Fourier based methods 
for estimating the mean and width, it can be shown that the 
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pulse-pair method is superior in cases where the 
signal - to -noise ratio, SNR, is small [6] . The pulse-pair 
method also offers an advantage in terms of calculation speed 
over Fourier based methods. 


F- Factor 

A hazard index or F- factor has been defined by Bowles 
[7] to measure the danger associated with microburst. The 
F- factor takes the form 

F = W' x /g - Wh/v (8) 

where w' x is the rate of change of the wind velocity along 
the flight path, v is the relative speed of the aircraft, g 
is the acceleration due to gravity, and is the vertical 
component of wind velocity along the flight path. A 
forward-looking Doppler radar can only measure radial 
velocities of wind along the flight path of the airplane and 
therefore estimates of can not be obtained. The radial 
component, Fr, can be written as 

F r = (v/g) (Aw x /AR) (9) 

where AW X is the change in radial velocity between range bins 
and AR is the distance between range bins. Estimates of the 
average radial veloctiy, W x can be obtained from the return 
data using the pulse-pair algorithm. 
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Statement of the Problem 

The ability of an airborne Doppler radar system to 
detect a microburst in the near terminal area depends heavily 
on the capabilities of that system to measure the horizontal 
component of windspeed versus range. The reflectivity 
associated with the precipitation present in the microburst 
is not the only contributing factor in this detection 
problem. Clutter return from objects in the near terminal 
area contributes to a loss in the signal - to - noise ratio and 
can impose a bias on estimates of windspeed. This thesis 
looks at reducing the effect of ground clutter on the 
pulse-pair mean estimates of windspeed through the use of 
clutter rejection filters based on modeling the clutter 
return in the near terminal area. 

The modeling techniques and coefficient estimation 
algorithms used in this thesis are discussed in Chapter II. 
Emphasis is placed on modeling the clutter as an 
autoregressive process and using an adaptive least squares 
prediction error lattice filter for determining the values of 
the model coefficients. Block processing and gradient 
adaptive algorithms are also discussed as alternative means 
for determining the model coefficients. 

In Chapter III, a model is used to generate microburst 
and clutter return data. The model was developed by NASA and 
uses actual clutter data taken from the Denver - Stapleton 
airport. Inputs to the model include microburst and radar 
system parameters. Using data generated from this model, 
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analysis of the pulse -pair estimates and power levels of a 
dry and a wet microburst are made in order to develop test 
cases for evaluating the performance of the adaptive clutter 
rejection filters. Spectral analysis and pulse-pair estimate 
analysis of the clutter returns are also given in Chapter 
III. The model clutter is observed to have the majority of 
its power contained in a narrow band of the Doppler spectrum 
and centered around zero. The clutter spectrum for this data 
set does not contain specular components as in the case of 
return from traffic on a nearby highway but this scenario 
should not be ruled out. 

The optimum clutter rejection filter is defined in 
Chapter IV, based on modeling the clutter in each range cell 
as an AR process. A complex normalized form of the least 
squares adaptive prediction error lattice filter, which 
offers possible fixed point implementation and involves fewer 
update equations than the unnormalized form, is used to 
determine the model coefficients. Analysis of pulse-pair 
estimates and power levels after filtering reflects the 
ability of the filters to remove the clutter without further 
biasing the pulse-pair estimates. Two other filtering 
schemes, one based on modeling the clutter to obtain a single 
filter for use over several range cells, and one based on 
conventional filtering techniques used in land based air 
traffic control radar systems, are also evaluated. 

Chapter V contains conclusions and recommendations for 
future work. Appendix A contains a derivation of a complex 
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least squares adaptive lattice structure. Appendix B 
contains a listing of the parameters used in the simulation 
model. Appendix C and D contain spectral estimates of dry 
microburst returns with and without clutter, respectively. 
Appendix E gives spectral estimates of the clutter return. 
Appendix F and G give the filter magnitude and phase 
response, respectively, based on modeling the clutter in each 
range cell. Appendix H, I, and J contain filtered spectral 
estimates where the filters were based on modeling the 
clutter in each range, modeling the clutter in range cell 20, 
and using a pulse canceller, respectively. 



CHAPTER II 


CLUTTER MODELING 
Introduction 

Data collected from a range cell may contain both 
weather return and ground clutter (e.g. return from ground 
structures) . Figure 6 represents a range cell containing a 
return from a microburst in addition to clutter. The clutter 
tends to bias the pulse -pair estimates of windspeed mean and 
width. Clutter may be centered around zero Doppler due to 
stationary objects on the ground or very specular in nature 
as in the case of traffic on a nearby highway. 

With clutter and weather spectra separated in frequency, 
a filtering scheme should be able to reduce the effects of 
clutter on the pulse-pair estimates of windspeed. Haykin [8] 
has suggested a means for modeling the clutter in the airport 
terminal area as an autoregressive, AR, process. The model 
can then be used to eliminate the clutter from the return 
signal. Gibson [9] has also used this modeling technique to 
design adaptive clutter rejection filters for use in an air 
traffic enviroment. The clutter rejection filters based on 
modeling the clutter offer a better estimate of the character 
of the clutter spectrum than the traditional moving - target 
indicator [9], MTI, filter. The MTI filter is a high pass 
filter in the form of a pulse canceller. This filter is used 




meters/sec 


Figure 6. Fourier spectral estimate of a range cell 
containing microburst and clutter return. 
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as a gross means for removing the clutter and does not take 
into account the spectrum of the clutter other than the fact 
that it is usually centered around zero Doppler. 

ARMA Modeling 

The AR process used for modeling clutter in the airport 
terminal area is actually an extension of the autoregressive - 
moving average, ARMA, model [10] used to describe a 
stochastic process. The ARMA model expresses the stochastic 
process in terms of a linear difference equation 


N M 

YT = S]<. Yp-k + £ d^ Ut-]c (10) 

k=l k=0 

where d(0) = 1. The input, u-p, to the system is assumed to 
be generated by a white noise process. The z- transform of 
the transfer function between the input, u-p, and the output, 
y-p, in Equation (10) can be written as 


H (z) 


D (z) 
A (z) 


(ID 


where 


N 

A(z) = 1 + Z a^ z'k (12) 

k=l 


and 


M 

D(z) = 1 + Z d^ z'^ . 
k=l 


( 13 ) 
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The power spectrum of the output of the ARMA model, P y (el w ), 
takes the form 

P y (e3W) = |H(e9W) |2p u(e jW) (14) 

where P u (e9 w ) and H(el w ) are z - transforms of the input and 
transfer function, respectively, evaluated at z = e9 w [11] . 

AR Modeling 

Wold's decomposition theorem [10] allows for an 
approximation of the ARMA model given in Equation (10) . The 
model of interest is an autoregressive, AR, model of the form 

N 

YT = - £ y T . k + u T (15) 

k=l 

with an associated transfer function 


H 


AR 


1 

A (z) 


(16) 


Wold's decomposition theorem states that an ARMA model of 
order (N,M) can be represented as an AR model of possibly 
infinite order. 

The power spectrum for the AR process in Equation (15) 
can be written as 


P 


AR 


(f) 


o 2 T 


N 

1 + Z a exp ( - j 27rfTk) 
k=l k 


2 


(17) 
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where a 2 is the variance of the input white noise process, 
and T is the sampling period. The input to the system was 
assumed to be generated by a white noise process to allow the 
frequency response of the input to be represented as a flat 
spectrum of magnitude a 2 . The power spectrum given in 
Equation (17) was also developed by Burg in maximizing the 
entropy associated with estimating unknown autocorrelation 
lags [12]. Park [13] discusses several advantages of Burg's 
maximum entropy method, MEM, in spectral estimation over 
Fourier based techniques. These advantages are revealed in 
its high resolution capability and absence of leakage. 


Levinson -Durbin Algorithm 

A method for determining the coefficients of the AR 
model can be obtained from the Yule-Walker equations 


N 

r y (k) = - E ajry(k-j) + cr 2 k = 0 
3=1 


N 

r y (k) = - E a-jr v (k- j ) 1 < k < N 

j=l 


(18) 


which relate the coefficients of the model to the 
autocorrelation of the process, r y (k) . Equation (18) can be 
written in matrix form as 
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r y (0) r y ( - 1) 
r y (1) r y (0) 

r y (N) r y (N- 1) 


r y ( -N) 
r y (- (N-l) ) 


r y (0) 





— — 


1 


02 


a l 

= 

• O 


a N 


0 


— — 




(19) 


Estimates of the negative autocorrelation lags can be 
obtained from the fact that r* y (m) = r y (-m), m = 0,...,N, for 
the stationary data case. The matrix on the far left of 
Equation (19) is defined as Toeplitz due to the equal- valued 
elements along each diagonal of this Hermitian matrix. 
Equation (19) containing the Toeplitz matrix can be solved 
using the Levinson -Durbin Algorithm [10] which provides a 
recursive method for determining the coefficients, {a -±, a. 2 , 

. . . . a^} • This algorithm is given in Table I. The 
Levinson-Durbin algorithm provides coefficients, [a^p o2p] 
... ( a N, 1 - • • a N, N , for all orders of the model up to N. 


Linear Prediction Error Filter 
The AR model in Equation (15) can also be interpreted as 

a linear prediction error filter, LPEF, that estimates y-p 
based on previous samples of the output with an associated 
estimation error, u T [11] . The output from the prediction 
error filter can be considered as samples from a white 
process due to the definition of the input, upi, in the AR 
model. Therefore, the filter in Figure 7 is sometimes termed 
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TABLE I 

Levinson-Durbin Algorithm 


o 2 o = r y {0) 

a l, 1 = -r y (l) /r y (0) 

for i = 1,2, ... ,N 
i - 1 

k i = - [r y (i) + l ai-i jr y (i- j) ] fo 2 x-l 

3=1 

a i , i = kj_ 

a i , j = a i-l, j + k i a *i-l ( j-l 1 < j < i-1 

a 2 i = (1 - |kj_ \ 2 )a 2 i . 1 

where 


a j “ a N , j 


1 < j < N 
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a whitening filter. The resulting transfer function for the 
whitening filter 

H w (z) = A (z) (20) 

is just the inverse of that for the AR model in Equation 
(16) . This whitening filter is a finite impulse response, 
FIR, filter which lends itself to the tap delay line 
implementation as seen in Figure 7. 

Lattice Prediction Error Filter 
The Levinson -Durbin algorithm also provides insight into 

another implementation of the linear prediction error filter 
defined in Figure 7. This new implementation involves two 
different errors, a forward linear prediction error, e N >p, 
and a backward linear prediction error, r^^, where 

N 

e N,T = £ a N,mYT-m (21) 

m=0 

N 

r N, T = £ b^my-T-m (22) 

m=0 

and ajj,o = bN,N = • The forward linear prediction error in 

Equation (21) is the error in predicting y-p based as N past 
samples, Cyt-I' Yt - 2 ' • • • ' Yt-N^- T b e backward linear 
prediction error in Equation (22) is the error in predicting 
Yt-N based on N future samples, {yt> Yt-1'---' Yt-N+1^- 

Haykin [11] shows that for a stationary process the 
coefficients for the backward prediction error are the 
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complex conjugates of the forward prediction error 
coefficients , 


b N,m “ a *N,N-m+l m = 1/2, ...,N (23) 
in reverse order. 

Based on the Levinson - Durbin algorthim, an order 
recursive implementation of Equations (21) and (22) can be 
formulated using the reflection coefficient, k n +i , [11] where 

e n+l , T = e n, T + k n+l r n,T-l (24) 
r n+l,T = r n,T-l + ^ n+l e n,T (25) 

for n = 0,...,N-1. This order recursive implementation of a 
linear prediction error filter can be implemented as a 
lattice filter, shown in Figure 8. 

The lattice prediction error filter has several 
properties which make it advantageous over the direct 
realization of a FIR filter. Haykin [11] shows that the 
reflection coefficients, k n+ ^, must have magnitude less than 
one and that this assures minimum phase for the filter (no 
zeros outside the unit circle). Also, Friedlander [14] 
refers to several studies which favor a lattice filter over a 
direct realization in terms of the effects of roundoff noise 
due to finite word length. Another major advantage is the 
ability to obtain the outputs from the lattice filter for all 
orders up to N using only N reflection coefficients. A 
direct realization would require N 2 -2N coefficients in order 
to implement the N required filter outputs. 



Coefficient Determination 
Burg ' s Block Implementation 
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The Burg algorithm [11] defines a block implementation 
for determining the reflection coefficients from the actual 
sampled data, y T . Once the reflection coefficients have been 
determined, the Levinson - Durbin algorithm can be used to 
compute the filter coefficients. Burg's method is based on 
minimizing the sum of the forward and backward prediction 
error energies. This cost function takes the form 

M-l M-l 

En - X I e n, T I 2 + l I r n, T I 2 (26) 
T=n T=n 

where M is the data length. This function is minimized with 
the constraint that the Levinson - Durbin algorithm holds. 
Taking the derivative of E n with respect to k^ , i = 
1,2,...,N, and setting the result equal to zero yields [15] 


k, 

1 


M-l * 

2 T | i C i-l,T-l e i-l,T 


M-l 

E 

T=i 


( Ir i-1,T-1 



|e i-l,T 


2 

I ) 


(27) 


The filter coefficients, for all orders up to N, can be 
obtained by recursively implementing the Burg algorithm in 
Equation (27), the Levinson -Durbin algorithm in Table I, and 
the forward and backward error order updates in Equations 
(24) and (25) . The block implementation of the Burg 
algorithm requires large amounts of storage space and 
considerable computation. 
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Adaptive Implementation 

An adaptive approach to determining the reflection 
coefficients reduces the large storage requirements needed to 
process the data as a block of information. Two methods 
which have found considerable use are the gradient adaptive 
approach and an adaptive scheme based on least squares. 
Besides reducing the storage requirements, the adaptive 
schemes are able to track changes in a nonstationary 
environment . 

Gradient Adaptive 

A complex gradient adaptive algorithm presented by 
Symons [16] is based on the cost function 

T 

H n+1 , T = £ M e n+l,il 2 + I r n+l , i I 2 J < 28 ) 

i=l 


where 


e n+l,T - e n, T * k n+l, T r n, T- 1 (29 > 

and 

r n+l , T = r n, T- 1 * k *n+l,T e n,T • ( 3 °) 

The gradient of the cost function with respect to k n+l,T 
yields 


5h 

6k 


n+1 , T 
n+1 , T 


2 I [e „ .r . , + r* .e 
i=l n+l,i n,i-l n+l,i n 



( 31 ) 
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An estimate of the gradient in Equation (31) can be obtained 
by evaluating it at (i = T) to form an instantaneous estimate 
[11] . The time update of the estimated reflection 
coefficients, k n+l,T' can now be written as 

k n+l,T+l = k n+l,T ' Mn, T [ e n+l , T r *n, T- 1 

+ r *n+l , T e n, (32) 

where Mn,T i- s tke step size controlling the convergence 
properties of the algorithm. Griffiths [14], [17] has chosen 
Mn,T to be inversely proportional to the prediction error 
power, R n ,T* Table II gives one version of the gradient 
adaptive implementation given by Giffiths which has been 
modified based on Symons results to support complex data. 
The associated lattice structure is equivalent to the one 
given in Figure 8, but with the reflection coefficient, k n +p, 
replaced with the reflection coefficient, * k n+l,T- 
Friedlander [14] provides a means for determining the AR 
model coefficients from the estimated reflection 
coefficients, shown in Table III, which is based on the 
transfer functions, Ap^-p^ 2 ) and ®N,t( 2 ^, between the input x<p 
and the errors, epj ( T and rp^T, respectively, where 


An, t 1 

(z) = 1 + a N/ pz" 1 

+ a N,2 2 ’ 2 + • • 

• • + a N,N 2 " N 

(33) 

B N, T 

(z) = 1 + b N> pz‘ 1 

+ b N ,2 2 ' 2 + • • 

• • + b N,N 2 " N • 

(34) 


This lattice structure implemenation of the transfer 
functions is given in Figure 9. 
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TABLE II 

Gradient Adaptive Lattice Algorithm 


Initialization : 

K p, -1 = 0 P = 1/ • • • /N 

r p,-l = r p, -2 = e p,-l = 0 P = N 

Rp f . 1 = 0,0 is an a priori estimate of the error 
covariance 

For each data point 

e 0,T = r 0 , T = YT 
Do for p = 0 N-l 

k p+l,T = k p+l,T-l + P t e p+l ,T-l r *p-l,T-2 
+ r *p+l,T-l e p,T-ll / R p,T-l 
e p+l , T = e p,T ‘ k p+l,T r p,T-l 
r p+l,T = r p,T-l ' k *p+l,T e p,T 

Rp , t = + f e2 p , T + r2 p,T-l3 

P: controls convergence rate 

a: reduces influence of past data samples 




31 


TABLE III 

Algorithm for Determining the AR Model Coefficients 
Based on a Gradient Adaptive Lattice Structure 


Initialize: = 0, p = 0,...,N-1 

For i = 0 , . . . ,N: 

ag ( i = bQ f i = [1 for (i = 0), 0 for (i > 0)] 
For p = 0 , . . . , N- 1 : 

a p+l,i = a p,i " kp+lbp,i-l 
bp+l,i = bp f j_-i - k p+iap,i 
a i = a N, i 
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Least Squares Adaptive 

An exact solution for the reflection coefficient 
order -time update can be found using a least squares, LS, 
adaptive implementation. In the LS case, the sum of the 
squared magnitude of the forward error, up to time T, 

T 

6 = Z I 1 1 2 (35) 

t=0 

is minimized with respect to the coefficients i = 

1,...,N} using the principle of orthogonality. The resulting 
LS adaptive lattice implementation of the prediction error 
filter [10] is given in Table IV and the associated lattice 
filter in Figure 10. A derivation of this LS adaptive 
lattice algorithm can be found in Appendix A. The AR model 
coefficients can be obtained in a similiar manner to that for 
the gradient adaptive algorithm. Friedlander [14] provides 
the algorithm in Table V and the associated lattice structure 
in Figure 11. A derivation of this lattice structure 
implementation of the transfer function is given by Honig and 
Messerschmi tt [18] for the real coefficient case. A 
derivation for the complex coefficient case is given in 
Appendix A. A major advantage of the LS adaptive 

implementation over gradient adaptive methods lies in the 
superior convergence rates of the LS implementation due to 
the exact solution it gives for the adaptation [14] . 
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TABLE IV 

Complex Least Squares Adaptive Lattice Algorithm 


Input parameters: 

N = maximum order of lattice 
Yt = data sample at time T 
a = exponential weighting factor 
Do for T = 0 to T max 
e 0,T = r 0,T = X T 

rS 0,T = Rr 0 , T = OR e 0,T-l + YtY*T 

yo,t = 1 

Do for n = 0 to {min(N,T)-l} 

k n+l,T = ^n+l , T - 1 + e *n, T r n, T - l/Yn, T- 1 

Yn+1 , T “ Y n , T * r *n, T R ' r n, T r n, T 

kr n+l,T = k n+l , T R " r n, T - 1 

e n+l,T = e n, T * k * r n+l , T r n, T- 1 

RS n+l,T = RS n, T ' kr n+l,T k *n+l,T 

ke n+l,T = R " e n, T k n+1 , T 

r n+l , T = r n, T- 1 * ke n+l,T e n,T 

Rr n+1,T = Rr n, T- 1 ' k *n+l,T ke n+l,T 

Note: Division by zero where y = y, R r , R e , set 1/y = 0. 

Initialize the variables k, r, R e , R r , and y to zero. 
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TABLE V 

Algorithm for Determining the AR Model Coefficeints 
Based on a Complex LS Adaptive Lattice Structure 


For i = 0,...,N do the following 
b p,-l,T = 0 f° r P = 0 , . . . ,N- 1 
a 0 , 0 , T = b 0,0,T = 1 for i = 0 
c 0 , 0 , T = 0 for i = 0 

a 0 , i , T = b 0 , i , T = c 0 , i , T = 0 for i > 0 
For p = 0, . . . ,N-1 

b p , i , T - 1 = b p , i , T + c p,i,T r p,T/Yp,T 
c p+l,i,T = c p, i , T + b p, i, T R ' r p,T r *p,T 
a p+l , i , T = a p, i , T ' b p, i-l,T-l R " r p,T-l k *p+l,T 
b p+l , i , T = b p, i - 1 , T- 1 ' a p, i,T R ' e p,T k p+l,T 
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Model Order 

In the AR modeling process the appropriate model order N 
is usually not known a priori. The use of too low a model 
order results in a smoothed spectral estimate and too high a 
model order induces spurious detail into the spectral 
estimate. A method for determining the correct model order 
based on some error criterion can be devised. The 
coefficient estimation techniques presented previously have 
monotonically decreasing error powers which prevent the use 
of the error power as the sole determining factor of the 
correct model order. Akaike [19], [20] has developed two 
criteria for determining model order selection. The first 
criterion is based on minimizing the average error, Ej^, for a 
one step prediction and is termed a final prediction error, 
FPE, criterion. The FPE criterion takes the form 


fpe n - ^ 


M+N+l ' 
M-N-l 


(36) 


where M is the number of data samples used. The model order, 
N, is choosen so as to minimize FPE. This method has been 
found to predict too low a model order in some cases. Akaike 
has developed a second criterion which is based on minimizing 
an information theoretic function assuming the process has 
Gaussian statistics. The Akaike information criterion, AIC, 
takes the form 


AIC n = ln(E N ) + 2(N+1)/M 


( 37 ) 


where the model order, N, is choosen to minimize the AIC. 



CHAPTER III 


THE EFFECTS OF CLUTTER IN 
ESTIMATING WINDSPEED 

NASA Model 

To support the development of an airborne pulsed Doppler 
radar system for the detection of microbursts in the near 
terminal area, NASA [2] has developed a simulation model for 
generating pulsed Doppler weather radar data containing 
microburst and clutter returns. The model is designed to 
allow for the input of clutter and microburst information as 
well as radar system parameters. Appendix B gives a listing 
of the possible input parameters to the model. 

To study pulse-pair estimates of radial wind velocity, 
two data sets were generated using the simulation model, one 
set contained a wet microburst and one set contained a dry 
microburst with the microburst centered 5km from touchdown. 
Appendix C contains Fourier spectral estimates of the 
simulated dry microburst and clutter return. Each data set 
included clutter based on actual clutter collected in the 
vicinity of the Denver Stapleton airport. In each case the 
data represented a stationary situation with the the aircraft 
positioned 7km from touchdown on a 3 degree glide slope with 
the first range bin recorded 1km from the aircraft. 
According to the simulated radar pulse width, successive 
range bins were located 150m from each other. Radar returns 
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were collected over 40 range bins from 512 pulses transmitted 
by the radar at a pulse repetition frequency of 3723.0 Hz. 
Other parameters related to the model for both the dry and 
wet microburst can be found in Appendix B. 

Bias in Pulse- Pair Mean Estimates 
Due to Clutter 

Figure 12 shows the resulting pulse-pair estimates of 
average windspeed for each range bin in the case of a radar 
return containing a wet microburst plus clutter. Equation 
(1) was used to convert the pulse-pair mean frequencies to 
the corresponding Doppler velocities. This figure 
illustrates the characteristic "S" curve associated with 
microburst. Estimates of average windspeed go from positive 
to negative as one moves through the microburst with zero 
radial velocity near the center of the microburst, range cell 
27. In the case of the dry microburst in Figure 13, the 
familiar "S" curve is no longer distinguishable. Two 
factors, bias due to clutter [21] and low signal - to -noise 
ratios, SNR, contribute to these poor estimates of windspeed 
within the dry microburst by the pulse-pair estimator. The 
effect of the bias due to clutter is less evident in the case 
of a wet microburst due to the expected higher 
signal - to -clutter ratios, SCR's. Due to the overwhelming 
effect of clutter and low SNR's on the pulse-pair estimates 
of windspeed, the dry microburst will be choosen as the test 
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case where any findings should also apply for the wet 
microburst case with higher SNR's and SCR's expected in a 
majority of the range cells. 

Analysis of a Dry Microburst 

A clutter rejection filter has been proposed as a means 
for removing the bias in the pulse-pair estimates. Before 
attempting to define this clutter rejection filter, a study 
of the pulse-pair estimates of a dry microburst without 
clutter is essential for interpreting the results after 
filtering. In the dry microburst case, the windspeeds versus 
distance from the of the center microburst are still 
characterized by the familiar "S" curve. Pulse-pair 
estimates in the presence of the dry microburst without 
clutter should therefore reflect this physical phenomena. 

The model is now used to generate the same dry 
microburst but without clutter. The pulse-pair mean 
estimates for this set of data are shown in Figure 14. Range 
cells 23 through 32 in Figure 14 are indicative of the "S" 
curve associated with the microburst. A spectral analysis, 
using a 512 point DFT, of range cells 25, 27, and 29, shown 
in Figure 15, reveals that these range cells do indeed 
contain microburst information with SNR's on the order of 20 
dB. As a note of reference, all Fourier spectral analysis 
presented in this thesis will be based on 512 data samples 
with no zero padding and will be normalized by the number of 
data points. 




meters/ s 



Figure 14. Pulse-pair mean estimates of the return from 
a dry microburst . 
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In Figure 14, range cells outside of 23 through 32 have 
what appear to be random estimates of mean windspeed. Figure 
16 provides a spectrum of range cell 15 and is typical of the 
spectral content in these range cells where the SNR's are 
very low. This low SNR can be used to explain the random 
mean estimates of windspeed found in these range cells [6] . 
Spectral estimates for all 40 range cells are given in 
Appendix D for the dry microburst without clutter case. 

Validation of Mean Estimates 
Width Estimates 

A method for determining valid pulse-pair mean estimates 
may be found by examining pulse -pair width estimates and 
power levels. Large widths estimates may indicate the 
presents of a bi -modal distribution or low SNR and serve as 
an indicator of questionable estimates on a range cell by 
range cell evaluation. This concept is supported by the 
width estimates of the dry microburst data without clutter 
given in Figure 17 . The large width estimates (approximately 
13 meters/sec) in the range cells away from the center of the 
microburst (range cell 27) are indicative of the low SNR 
expected in the dry microburst environment. 

Power Levels 

Power levels may also provide a level of significance 
that can be associated with the pulse-pair mean estimates of 
windspeed. Figure 18 gives the power levels associated with 
the dry microburst data without clutter. The power levels 
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were computed as a sum of the squared magnitudes of the IQ 
radar return samples. The resulting power levels in this 
thesis are therefore normalized by the sampling period. The 
range cells near the center of the microburst (24-30) have 
approximately a lOdB to 20dB relative gain in power compared 
to range cells away from the center of the microburst. This 
gain can be associated with the higher target concentrations 
near the center of the dry microburst. The low levels of 
power indicate the presence of noise only. 

Threshold Levels 

The analysis given in the previous two sections suggests 
that the power levels and the pulse pair width estimates may 
be useful in forming a statistic to identify the levels of 
significance to associate with the pulse-pair mean estimates. 
The threshold power levels can be obtained based on 
statistics of noise and clutter power levels. After 
filtering, the power in the return data can be contributed 
two sources: noise and weather information. Thresholds for 
width estimates may be defined by statistical analysis of 
radar returns in windshear situations. The actual statistic 
will not be formulated in this thesis, but the power levels 
and pulse-pair width estimates will be used as guides for 
evaluating the significance of pulse-pair mean estimates. 
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Clutter Analysis 

At this point, a study of the clutter may reveal 
possible filtering techniques for removing the bias imposed 
by the clutter on the pulse-pair estimates of windspeed. The 
radar simulation model can be used to generate range bins 
containing clutter only. Pulse-pair estimates of mean 
clutter speed, where "clutter speed" refers to the Doppler 
spectrum associated with the clutter, are given in Figure 19. 
The clutter tends to be centered around zero Doppler due to 
stationary objects along the flight path, but a slight shift 
in speed to negative Doppler occurs in the first few range 
cells. This shift can be contributed to returns from 
sidelobes which contain information from objects at close 
range which have a different relative velocity with respect 
to the aircraft than objects along the intended flight path. 
The demodulation of the return signal to Doppler frequencies 
is based on the relative ground velocity of the aircraft 
along the intended flight path. 

Figure 20, containing range cells 6, 20, and 35, is 
typical of the spectral estimates of the clutter for this set 
of 40 range cells. Appendix E contains Fourier spectral 
estimates of the clutter for all 40 range cells. This 
particular set of clutter return data does not contain any 
specular clutter away from zero Doppler which may result from 
traffic on a nearby highway, but this scenario should not be 
ruled out. Even though the clutter tends to be centered 
around zero Doppler, the actual width and power associated 
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with the clutter varies from one range cell to the next. 
This variation can be seen in Figures 21 and 22 which 
represent the power and width estimates, respectively, 
associated with each range cell. These width estimates were 
obtained using the pulse-pair algorithm. The large width 
estimates in some of the range cells may be the result of low 
power levels within the these range cells which effect the 
performance of the pulse-pair estimator [6] . The actual 
clutter will also vary from airport to airport and even from 
runway to runway due to the varied physical environments. 
This imposes a problem in designing an optimum filter to 
remove the clutter from the radar return data. 
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Figure 21. Pulse-pair width estimates of the return 
from clutter. 



Figure 22. Power levels associated with the return 
from clutter. 





CHAPTER IV 


CLUTTER REJECTION FILTERS 

The Optimum Clutter Rejection Filter 
Based on observations of the variation in the simulation 
model clutter spectra over many range cells it appears that 
the optimum filtering scheme for removing the clutter from 
the return data in a particular range cell would need to be 
based on the geographical location of the range cell. Even 
for a given runway, periodic updates of the filters used on 
final approach would need to be obtained due to glideslope 
flight path variations and also changes in the physical 
environment of the airport terminal area. This scenario 
suggests some type of adaptive modeling process to define the 
clutter in each range cell with updating as warranted. 

Adaptive modeling can be accomplished using the AR model 
defined in Chapter II implemented as an adaptive LS lattice 
structure in order to determine the model parameters. The 
lattice structure was chosen for its ability to provide 
coefficients of the model recursively for all orders up to N. 
This provides the ability to use different model orders for 
each range cell based on the requirements mandated by the 
clutter. The adaptive LS lattice structure implementation 
for determining the model coefficients is choosen for its 
superior convergence rates over gradient adaptive schemes. 
The lattice structure also insures that the resulting filter 
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will be stable. The model parameters could be obtained on 
approach to a runway on days when favorable weather 
conditions exist. In the presence of a microburst or other 
source of windshear, the best possible filters could then be 
available to remove unwanted ground clutter from the return 
data, enhancing the ability to detect windshear. 

The Complex Square Root Normalized Recursive 
Least Squares Lattice Estimation Algorithm 

An alternative form of the adaptive LS lattice structure 
given in Chapter II can be obtained, involving fewer update 
equations and potential for fixed point implementation. This 
alternative implementation is termed a normalized recursive 
LS lattice structure. A complex form of the algorithm, 
needed to process the I and Q data, was not available in the 
literature. A modification of the normalized recursive LS 
lattice algorithm for real data given by Lee [22] has been 
developed to handle complex data. The derivation is given in 
Appendix A. Table VI gives the necessary equations needed to 
implement the complex square root normalized recursive least 
squares lattice estimation algorithm given in Figure 23. 

The complex square root normalized recursive least 
squares lattice estimation algorithm requires only 3 
equations per order update compared to six in the 
unnormalized version in Table IV. The normalizations also 
restrict the magnitude of the normalized reflection 
coefficients, T, forward errors, v, and backward errors, r| to 
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TABLE VI 

Complex Square Root Normalized Recursive Least 
Squares Lattice Estimation Algorithm 


Initialize : 

r_-L = o a is a small positive value 

for T = 0,T max 
R t = cxRt- 1 + y*TVT 

v 0 , T = no,T = Yt( r T^ * 1//2 

for n = 0 to [min{T,N} -1] 

r n+l,T = d ' v *n, T v n,TJ 1//2 r n+l,T-l 

C 1 ‘ T l*n,T-l T ln / T-l] 1//2 + v *n,Tn n ,T-l 
v n+l , T = t 1 * r *n+l,T r n+l,Tl ‘ 1/f2 [v n>T - r * n+1 , T n n , T- d 
d • T l*n, T- l’ln, T- 1^ ‘ 1 ^ 2 

*ln+l , T = d - r* n+ 1 /T r n+1;T ] - 1 / 2 [Tin, T- 1 • r n+l,T v n,T^ 
d - v* n(T v nfT ] - 1 / 2 

Note: Division by zero where y = 1/x : x = 0 should result in 
y = 0. Initialize the variables T, v, and r| to zero. 
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values less than one. This conditioning of the variables 
allows for fixed implementation. 

The transfer function for the normalized recursive least 
squares lattice estimation algorithm in Figure 24 can be 
written as 


A N( t ( z) 

= ^N, 0 

B n , t ( z) 

= &N,0 

C N , T (z) 

= SN, 0 


l 2 " 1 + Sn, 2 z ' 2 + 
32n , 1 z ~ 1 + , 2 2 2 + 

£N , 1 z " 1 + £N,2 Z ' 2 + 


• • + a.N, N z (38) 
. . + N z ( 39 ) 

• • + c N(N z* n . (40) 


The normalized transfer functions An,T(z) and , T ( 2 ) 
correspond to the transfer functions between the normalized 
input, y-pR-p' 1 / 2 , and the normalized forward, v N/Tj and 
backward, t)n,T, errors, respectively. The normalized transfer 
function Cj^p-tz) corresponds to an auxiliary variable 
defined in Appendix A and is only needed for the 
implementation of the normalized transfer functions Apj^>p(z) 
and B N , T ( Z ) as a lattice structure. 

The coefficients for the transfer function of the AR 
model defined in Equation (15) in Chapter II, and equally so, 
those for the transfer function of the whitening filter 
defined in Equation (20) of Chapter II, can be obtained from 
the normalized transfer function given in Equation (38) as 


= (aN,0/^N,0) + (^N,1/^N,0’ 
+ ... + (a NfN /a N(0 )z' N 


\ j\j , z / —iN , 


:'l + a™ oz - 2 + 


An algorithm for obtaining the normalized coefficients is 


INPUT 


SECTION 

1 


—1 T ^ z ^ 


En- 1,T (z) section 

>- N 

%-l,T (z) 


(a) The overall structure 


A single section. 


Figure 24. Transfer function implemented as a lattice 
structure for the complex square root normalized least 
squares adaptive prediction error filter case. 
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given in Table VII. This algorithm is given by Honig [18] 
and Fr iedlander [14] for the real data case. A derivation 
for the complex case is given in Appendix A. 


10th Order AR Clutter Model 

Spectral estimates of the range cell clutter data using 
the power spectrum of the AR process in Equation (17) can now 
be compared to Fourier estimates to evaluate the ability of 
the AR coefficients to model the clutter. Based on 
preliminary evaluations a model order of 10 was chosen for 
the AR process in modeling the clutter. This fixed model 
order simplifies the modeling process with the constraint 
that the model order is sufficient in a majority of the range 
cells. Future evaluations based on Aikaike's criterion in 
Chapter I or other model order determination techniques may 
be used if the exact model order is vital for each range 


cell . 

Figure 25 contains a Fourier spectral estimate, based on 
a 512 point DFT , of the clutter contained in range cell 6. 
The ability to model the clutter with a 10th order AR model 
is evident by comparing the model given in Figure 26 with the 
corresponding Fourier spectral estimate. Note the ability 
to model the clutter in Figure 26 even when the mode of the 
clutter is shifted from zero Doppler. This ability to model 
the clutter as a low order AR process satisfies the necessary 
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TABLE VII 

Algorithm for Determining the AR Model Coefficients Based 
on a Normalized Complex LS Adaptive Lattice Structure 


R-p = ocRt- l + 

For i = 0 , . . . , N do the following 
bp , - 1 , T = 0 for p = 0 , . . . , N- 1 

a 0,0,T = b 0,0,T = %' 1/2 for 1 = 0 
c 0 , 0 , T = 0 for i = 0 

ao, i,T = b 0,i,T = c 0,i,T = 0 for 1 > 0 

For p = 0, N-l 

b p , i , T - 1 = C b p , i , T + c p ,i, T n p ,T] 

[1 ‘ T l*p, T^p, * 1//2 

Cp+i,i, T = C c p, i , T + b p,i,TH*p,T] [l-n*p,Tnp,Tl 1/2 
a p+l, i,T = t a p,i,T ' b p,i-l,T-l r P+1 ,t1 

[i-r*p+i, T r p + 1/T ] ' 1/2 

b p+ 1 , i , T = [ b p , i - 1 , T - 1 * a p, i,T r p+l,T^ 

u-r* p+ i,T r p+i,T] ‘ 1/2 

where 

a p , i = a p,i/ a p,0 
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limitations on storage requirements in a real world 
implementation, in that only a limited number of coefficients 
need be stored and completely specify the model. 

IQth Order Complex Coefficient FIR Filters 
The coefficients of the AR model can now be used to 

filter the microburst plus clutter data. The whitening 
filter defined in Equation (20) of Chapter II gives the 
relationship between the model parameters and the filter 
coefficients and can be used as a clutter rejection filter. 
The discrete time filter takes the form 

x T = y T + a N; iy T -i + ... + a N;N y T . N (42) 

where y T is the radar return containing microburst plus 
clutter data and x>p is the result of the filtering process. 
This filter is a finite impulse respone, FIR, filter which 
can be implemented as a tapped delay line as seen in Figure 
27. Another advantage of this type of filter is the complex 
coefficients obtained from the modeling process, which define 
the filter unit sample response as complex, provide the 
ability to define filters with nonsymmetric magnitude 
frequency responses. This is necessary when the clutter is 
shifted away from zero Doppler or in the case of specular 
clutter which may be due to moving objects on the ground. 
The phase associated with this type of complex coefficient 
FIR filter is not be constrained to linear. A linear phase, 
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however, is not necessary for this filter if the pulse pair 
algorithm is used to estimate windspeed mean and width in 
each range cell [21] . 

Analysis of Filtering Schemes 

Filters Based on the Clutter 
in Each Range Cell 

Using the filter defined in Equation (42) , an analysis 
will be made of the effects of filtering each range cell with 
a 10th order complex coefficient FIR filter where the 
coefficients are determined from the clutter using the 
complex square root normalized adaptive least squares lattice 
estimation algorithm. This analysis will be based on 
simulated radar returns from a dry microburst plus clutter 
data where the effect of bias on the pulse-pair mean 
estimates of windspeed is significant. Any resulting 
generalizations from this analysis can be applied in the wet 
microburst case where the effect of bias is reduced due to 
higher SNR's and SCR's. The magnitude and phase response of 
the filters defined for each range cell are given in Appendix 
F and G, respectively. The spectral content of the filtered 
range cells is given in Appendix H. Figure 28 shows the 
pulse-pair estimates of the post filtered data containing a 
dry microburst plus clutter. The characteristic "S" curve is 
now evident in range cells 23 through 32. The ability to 
remove the clutter without further biasing the pulse-pair 
estimates of the dry microburst is evident in Figure 29 by 
comparing the post filtered estimates of windspeed to those 




range cell number 


Figure 28. Pulse-pair mean estimates of filtered dry 
microburst plus clutter data where the appropriate, filter 
based on the clutter in each range cell has been applied. 



□ — 

filtered estimates 

— 

dry microburst estimates 


Figure 29. Pulse-pair mean estimates of filtered dry 
microburst plus clutter data where the appropriate filter 
based on the clutter in each range cell has been, applied 
compared to estimates given in the case of a dry microburst 
without clutter. 
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from the data containing a dry microburst only. The 
fluctuating estimates of windspeed outside of range cells 23 
through 32 in Figure 28 may be attributed to instabilities in 
the pulse-pair estimator when the resulting filtered data 
contains a very low SNR. A low SNR after filtering is 
illustrated in Figure 30 with the filtered spectrum of range 
cell 7 . 

The pulse-pair mean estimates after filtering will now 
be evaluated using pulse-pair width estimates and power 
levels. Using the pulse-pair width estimates of the dry 
microburst without clutter, shown in Figure 17, as a basis 
for setting thresholds on valid mean estimates, assume that a 
width estimate threshold of 7 m/sec is set. This choice of 
threshold is not based on a statistical evaluation but is 
used only for relative comparisons. Pulse-pair width 
estimates from filtered dry microburst plus clutter data, 
illustrated in Figure 31, reveal that range cells 23-26 and 
28-31, which contain valid mean estimates, have width 
estimates below this threshold. This correlates with the 
previous estimates of returns from the dry microburst without 
clutter. With this criterian, the mean estimate at the 
center of the microburst (range cell 27) would be considered 
questionable. A large width estimate of approximately 
llm/sec is likely since the clutter and microburst 
information would have the highest probability of occupying 
the same frequency ranges near zero Doppler, and would both 
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be extracted in the filtering process resulting in a very low 
signal to noise ratio. 

A minimum threshold for power level estimates might also 
be defined using the dry microburst without clutter estimates 
in Figure 18. For comparison, a -105 dB threshold will be 
set. Based on this threshold, the power levels in Figure 32 
define range cells 23-31 and range cell 40 as having adequate 
return power levels to provide valid mean estimates. The 
power levels for range cells 27 and 40 do not correlate with 
the information provided by the width estimate analysis in 
labeling range cells 27 and 40 as having questionable mean 
estimates. This contradiction between width estimates and 
power levels, concerning valid mean estimates, reflects a 
need for comparing the results of one test against the 
results of another before making a decision concerning the 
validity of a mean estimate. 

Filter Based on the Clutter 
in a Single Range Cell 

The 10th order complex coefficient filters designed for 
each range cell were able to eliminate bias in the pulse -pair 
mean estimates due to clutter while preserving the microburst 
information. This ability to preserve the microburst 
information after filtering is evidence by the resulting 
output power levels and reasonable width estimates. An 
alternative approach to designing filters for each individual 
range cell is to define a single set of filter coefficients 
that can be used for several adjacent range cells or over the 



Figure 32. Power levels after filtering using the filter 
coefficients obtained from the clutter within each range 
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entire set. The design of this filter could still be based 
on clutter statistics. The filter designed for range cell 20 
was choosen to test the effects of using one filter over 
several range cells. The range cells selected for evaluation 
are those around the center of the microburst which contain 
clutter as well as microburst information. The magnitude and 
phase response of this filter is given in Figures 33 and 34, 
respectively. This filter has a stop band centered at zero 
Doppler which is characteristic of the clutter in the 
majority of these range cells as can be seen in Appendix F. 

Figure 35 gives the pulse-pair mean estimates after 
filtering with the clutter rejection filter designed for 
range cell 20. A comparison in Figure 36 with the mean 
estimates from the dry microburst without clutter seems to 
support the possible use of one filter over range cells 23 
through 32. The spectral content of the filtered range cells 
is given in Appendix I. In order to evaluate the mean 
estimates after filtering, an analysis of the power levels 
and width estimates will be made based on the thresholds set 
in the previous section. 

Figure 37 gives the pulse-pair width estimates for the 
filtered dry microburst plus clutter data. Based on a 
threshold of 7m/sec, range cells 23-26 and 28-31 would be 
defined as having valid mean estimates. This result compares 
well with those found using the filters designed for each 
range cell. The power levels resulting from filtering are 
given in Figure 38 and are consistent with those found in the 
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range cell number 


Figure 35. Pulse-pair mean estimates after filtering using 
the filter coefficients for range cell 20. 



Figure 36. Pulse-pair mean estimates after filtering using 
the filter coefficients for range cell 20 compared to 
estimates given in the case of a dry microburst without 
clutter . 





Figure 38. Power lev 
coefficients for range 
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previous case using the appropriate filter for each range 
cell. These results suggest the possible use of one filter 
over several range cells where the range cells are known to 
have a similar clutter spectrum. 

Comparison with a Pulse Canceller 

In order to evaluate the filters based on modeling the 

clutter, a comparison will be made using a pulse canceller 
filter which is found in land based air traffic control radar 
systems [9] . The pulse canceller takes the form of a first 
order difference equation with input y>p and output x<p where 

x T = y<p * Yt - 1 • (43) 

The magnitude and phase responses of this filter are given in 
Figures 39 and 40, respectively. Figure 41 gives the 
pulse-pair mean estimates after filtering with the pulse 
canceller. A comparison in Figure 42 with the mean estimates 
from the dry microburst without clutter case reveals a larger 
error in the mean estimates for range cells 23 through 31 
when compared to that found in the case of using the filter 
designed from the clutter within each range cell, Figure 28. 
The spectral content of the filtered range cells is given in 
Appendix J. For comparison, Figure 43 gives the filtered 
results for range cell 24 using the filter designed for range 
cell 24, the filter designed for range cell 20, and the pulse 
canceller, respectively. This comparison reveals a larger 
attenuation of the microburst information in the case of the 
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meters/ sec 

Figure 39. Magnitude response of a pulse canceller filter. 



meters/sec 

Figure 40. Phase respone of a pulse canceller filter. 
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Figure 41. Pulse-pair mean estimates after filtering using a 
pulse canceller. 




dry microburst estimates 


filtered estimates 


Figure 42. Pulse-pair mean estimates after filtering using a 
pulse canceller compared to estimates given in the case of a 
dry microburst without clutter. 



Rancre Cell 24 




meters/sec 

Filtered using range cell 24 model coefficients 


Rancre Cell 24 




meters/sec 

Filtered using range cell 20 model coefficients 


Range Cell 24 




meters/sec 

(c) Filtered using a pulse canceller. 

Figure 43. Spectrum of range cell 24 after filtering 
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pulse canceller. This is reflected in the lower power levels 
given in Figure 44 for the pulse canceller case near the 
center of the inicroburst. These resulting lower SNR's 
contribute to the large width estimates found after using the 
pulse canceller, shown in Figure 45. Only two range cells, 
28 and 29, have width estimates within the 7m/sec threshold 
previously set to define a valid mean estimate. 





CHAPTER V 


CONCLUSIONS 

The potential for disaster associated with low altitude 
windshear in the near terminal area mandates the need for 
some type of detection system. It is anticipated that an 
airborne pulsed Doppler radar system used to measure radial 
components of windspeed can help assess the hazard associated 
with windshear particularly in the case of a microburst. 
Clutter presents a major obstacle to the use of airborne 
radar at low altitudes as a remote sensor of windshear. 
Large clutter returns from objects in the terminal area can 
impose a bias on estimates of the radial components of 
windspeed. This thesis has addressed the use of adaptive 
clutter rejection filters to eliminate the source of this 
bias . 

The adaptive clutter rejection filters investigated here 
are based on modeling the clutter as a low order 
autoregressive process. A low order filter offers 
implementation advantages in terms of memory requirements and 
computational load, both advantageous for real time 
implementation. The adaptive property of these clutter 
rejection filters offers the ability to continuously update 
and improve the filters as new clutter return becomes 
available. Because the filter coefficients are obtained from 
a complex square root normalized recursive least squares 
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lattice estimation algorithm, stable clutter rejection 
filters are assured. The least squares adaptive method for 
determining the filter coefficients also offers a faster 
convergence rate compared to gradient adaptive methods, and 
the normalization allows for implementation using fixed point 
arithmetic further enhancing the potential for real time 
implementation. Fast convergence rates are important since 
the clutter return is generally non - stationary over the 
long-term but is assumed stationary over short observation 
intervals. Additionally, complex coefficient FIR clutter 
rejection filters allow for non - symme t r ic magnitude 
responses, required in cases where the clutter spectrum is 
non- symmetric about zero Doppler. 

Several filtering schemes used to eliminate the clutter 
have been considered. One such filtering scheme defines a 
model for each range cell at any point along the final 
approach to a runway. These models can then be indexed by 
geographical location for future use. It has been 
demonstrated here that the same model may be used for several 
range cells, thereby decreasing the amount of information to 
be stored. More research is needed to make general 
conclusions . 

The pulse-pair algorithm has been used herein to 
determine the average windspeed within a range cell from the 
filtered return data. The resulting mean estimates tend to 
be unstable or biased when the signal to noise ratio is 
small. This thesis has proposed the use of pulse-pair width 
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estimates and power levels to form a statistic for 
establishing confidence in pulse-pair mean estimates. 

Further research in the area of stationarity and 
repeatability as applied to the environment being modeled is 
warranted. Another area of research not addressed in this 
thesis but needing further study is the determination of the 
optimum model order. Further study of clutter environments 
is also needed before defining the appropriate filtering 
scheme to use. Particular emphasis will need to be placed on 
dynamic clutter models. 
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Appendix A 


The Complex Square Root Normalized Recursive 
Least Squares Lattice Estimation Algorithm 

Introduction 

The square root normalized recursive least squares 
lattice estimation algorithm presented by Lee, Morf, and 
Friedlander [22] was derived for the real data case. This 
appendix looks at adapting the algorithm to process complex 
data. The need for a complex form of the algorithm arises 
when processing IQ data from a radar return or other sources 
of data requiring quadrature sampling [16] . This appendix 
will not discuss the advantages of the square root normalized 
form of the algorithm, but references concerning the 
advantages are given by Lee [22] . Honig [18] and Friedlander 
[14] address the issue of obtaining the linear prediction 
error filter coefficients from the normalized lattice 
parameters for the real data case. A derivation will be 
developed for the complex coefficient case based on that 
presented by Honig. 

This appendix is divided into six major sections. The 
first section discusses the least squares approach to 
estimating the complex coefficients in an autoregressive 
model. The next two sections develop the order- and 
time-updates required in a recursive lattice estimation 
algorithm. A normalization is introduced in section four 
that reduces the number of updates required in the lattice 
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estimation algorithm. The fifth section is a summary of the 
equations and intializations required to implement the 
complex square root normalized recursive least squares 
lattice estimation algorithm. The final section uses the 
normalized lattice parameters to determine the complex 
prediction error filter coefficients. 


Autoregressive Model 
Least Squares Estimate 

The task of fitting an Nth-order autoregressive model to 
a complex set of data {x*-, 0<t<T} involves finding a complex 
set of coefficients {ajj A> to minimize the sum of 

the squared prediction errors 

I e I 2 = (x'* + A N '*X'*) (x + XA n ) (A . 1 ) 


where 


e = + XA n 

X' = [XQ, • • • ,X T ] 

An' = t a N, l - • • • ' 9 n,n1 

and 


0 

0 

0 

x 0 

0 

0 

- 

x o 

0 

• 

• • « * 

x 0 

- X T - 1 

X T - 2 • • • 

X T-N - 


(A. 2) 
(A. 3) 
(A. 4) 


(A. 5) 
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The transpose of a vector is noted as x’ / and the complex 
conjugate of a vector is noted as x* • The minimization is 
obtained by taking the partial derivative with respect to the 
complex coefficient vector A^ and setting the result equal to 
zero. The steps required for complex differentiation are 
discussed by Miller [23] . Expanding the product in Equation 
(A. 1) 

|e|2 = (x' *x) + (x' *XA^) + (An'*X'*x) + 

(A N '*X'*XA N ) (A. 6) 

and taking the partial derivative with respect to Ajj yields 
the following 


0 [x' *x] /3 Ajj = 0 

(A. 7) 

8[x'*XA N ] /3A N = X'x* 

(A. 8) 

9[A N '*X'*x] /3 A n = 0 

(A. 9) 

9CA N '*X'*XA N ]/9A N = X'X*A N * . 

(A. 10) 


Setting the sum of the partial derivatives in Equations 
(A. 7), (A. 8), (A. 9), and (A. 10) equal to zero and solving for 

A^ yields 

A N = - [X ' *X] -lx'*x (A. 11) 

the least squares estimate of the complex predictor 
coefficients . 


C 
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Projection Operators 

The complex coefficient vector in the prediction error 
in Equation (A. 2) can now be replaced by the least squares 
estimate of the complex coefficient vector in Equation (A. 11) 
to yield 

e = x - X [X' *X] ^X' *x . (A. 12) 

Using the least squares estimate of the prediction error in 
Equation (A. 12), a projection operator can now be defined as 

P N ■ X [X' *X] 'lx'* . (A. 13) 

This operator projects the vector x onto the subspace of past 
observations (the columns of X). Its orthogonal complement 
can be defined as 

PN 1 = (I - P N ) . (A. 14) 

The orthogonal complement of the projection operator P^ in 
Equation (A. 14) can now be used to define the prediction 
error 

e = P-*"nx (A. 15) 

as a projection of the observed data onto a subspace which is 
perpendicular to the one containing past observations. 

A property of projection operators that will be useful 
is P]sjPn = P N and = p_ *"N- A proof of this property 

of projection operators is 
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P N P N = X [X' *X] _1 [X' *X] [X'*X] -ix' * 

= X [X' *X] _1 X' * 

= P JyJ . ( A . 1 6 ) 

Order -Update Recursions 
Sample Space 

This derivation of the square root normalized recursive 
least squares lattice estimation algorithm will consider a 
single channel complex data case where Lee [22] applied it to 
a multichannel case for real data. An observation of time 
samples {x{- G C, 0<t<T} in vector form is defined as 

|x> T = [xq, xi, x T ] ' . (A. 17) 

Also, a delay operator, z‘ n , on |x>-p is defined as 

|z'1x>t * [0, xq , -•-/ x^ - i ] ' (A. 18) 

where x t , t<0 is assumed to be equal to zero. 

The linear space containing |x>-p is spanned by the 
T+l observations vectors { | x>-p/ I z " . . . , i z" T x>-p} • The 

vector inner product on PPp is defined as 

<x | y>T * ( |x>t' *) ( ly>T^ 

T 

= £ x t *y t (a. 19) 

t=0 


where |x>t, ly><p G . 

For convenience, the transpose of |x>>p is defined as 
T <x| . A matrix composed of elements of H<p will be defined as 
IX>t where the matrix transpose is noted as -j><x| > an< ^ t ^ ie 
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matrix inner product is defined as 

<X|Y> t ■ (|x> T '*> (|y> t ) . (A. 20) 

Subspace of Past Observations 

The subspace of H-p containing past observations is 

denoted by X^ nj T and is spanned by 

{ I z * ^"X> ’j' , |z' n x>i>)/ n < T (A. 21) 

where. 


l x l,n^T = [|z'^x><r, . . . , | z ' n x>-p] • (A. 22) 

The projection operator on the subspace Xp (11/ p is 

p l,n,T * l x l,n>T <x l,nl x l,n > T‘ 1<x l,nlT* • (A. 23) 

Also, the orthogonal complement of the projection operator 
defined on Xp^j-^p- is 

pJ -l,n,T = (I - Pl,n,T) • (A. 24) 

Coordinate Map 

In order to extract the most recent time sample, a 
coordinate map is defined such that n(|x>«p) = xp>. A vector 
form of this operator can be written as 

|n> T = [0,0, . . . , 0, 1] ' (A. 25) 

where ( T <n|) (|x>-p) = x<p . A coordinate projector can also be 
defined in terms of |n> T such that 
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P'p | x> ■j 1 - [0,0, . . . , x^] 


(A. 26) 


where 


P>p * | n> , p < Cn | . (A. 27) 

Forward and Backward Prediction Errors 

Projecting |x>'p onto the subspace of past observations, 

x l,n,T' yields an nth order forward prediction error 

I e n >T = ! X >T ' R l,n,Tl x ^T ~ P^"i , n , T I X ^T (A. 28) 

which is orthogonal to ; but lies in Xg ( n< t • The 

forward prediction error covariance can be defined as 

R e n ,T = <e n I e n^T = <e n |x>'p = <x|e n >>p . (A. 29) 

An nth-order backward prediction error vector can also be 
defined as 


|r n > T = lz’ n x> T - Po , n- 1 , T I z " n x>T 

= P^*0,n-l,Tl z " nx> T • (A. 30) 

A delayed version of the backward prediction error 

|z _1 r n > T = P- L l,n,Tl z ‘ n ‘ lx >T (A. 31) 

will be needed in the lattice recursions. Its covariance is 

Rr n,T-l = < z " lr n I z * lr n > T = <r n lr n >T-l • (A. 32) 

Decomposition of Subspaces 

The fact that |e n >-p lies in Xo /R/ t and is orthogonal to 
x l,n,T loads to the direct sum of the subspaces such that 
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x O,n,T = x l,n,T ® I e n>T • 


(A. 33) 


The projection operators can be updated in a similar fashion 
wi th 

p 0,n,T = p l,n,T + I e n > T <e n I e n>T 1<e n lT (A. 34) 

and its orthogonal complement 


p_ ^0,n,T = p_ *"l,n,T ' I e n > T <e n I e n > T 1<e nlT ■ (A. 35) 

Using the delayed backward prediction error, the projection 
operator, Pl,n,T' can updated as 

P l,n+1,T = p l , n , T + I z ' 1 r n > T < z ' 1 r n I z ‘ 1 r n > T ‘ 1 < z ' p r n I T * 

(A. 36) 


and its orthogonal complement as 

pi -l,n+l,T = P-h-^T - I z 1 r n > T < z ’^r n I z" 1 r n > T ' 1 <z" 1 r n I T *. 

(A. 37) 

Error Order -Update 

Equation (A. 37) is an expression for the projection 
order -update at time T. The order -update of the forward 
prediction error is obtained by operating on Ix>t using 
Equation (A. 37) to yield 

l e n+l>T = I e n>T ' I z ' 1 r n > T <z ‘ 1 r n I z ' 1 r n > T ' 1 <z ‘ 1 r n I e n > T . 

(A. 38) 


At this point, a reflection coefficient, k n +i >T , will be 
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defined as 

k n+l,T * < e n I z " lr n>T • (A. 39) 

Using the forward prediction error covariance in Equation 
(A. 32) and the reflection coefficient in Equation (A. 39), the 
forward prediction error order-update in Equation (A. 38) can 
be written as 

e n+l , T = e n, T ' r n, T- 1 R ' r n, T- l k *n+l , T (A. 40) 

while looking only at the Tth component. Comparing this to 
Lee [22] , a conjugate of the reflection coefficient in 
Equation (A. 40) is required for complex data due to the 
definition of the inner product in Equation (A. 19) and the 
reflection coefficient in Equation (A. 39). The covariance 
order -update for the forward prediction error in Equation 
(A. 38) is given by 

<e n+l l e n+l>T = <e n I e n>T * < e nl z ' lr n>T 

<z’ 1 r n | z' 1 r n > T ’ 1 <z' 1 r n |e n > T (A. 41) 

or 


Re n+1,T “ Re n,T ' k *n+l , T R ’ r n, T- l k n+l , T • (A. 42) 

Using Equation (A. 35), a similar order-update can be 
obtained for the backward prediction error. Operating on 
|z' n *lx>rp yields 

I r n+l>T = l z ‘ lr n > T ' I e n>T <e n I e n>T* 1 <e n I z ' lr n>T 


(A. 43) 
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or 

r n+l , T = r n, T- 1 ' e n,T R ’ e n,T k n+l,T • 

The covariance order -update for the backward prediction 
error is 

<^n+l I r n+l > T = <z‘ lr nl z " lr n > T ' <z ’ lr nl e n > T 
<e n l e n>T' 1<e nl z ' lr n>T 


or 


R r n+1 , T = Rr n, T- 1 - ^n+l , T* R ' e n, T^n+1 , T • (A * 46) 

Time -Update Recursions 
Decomposition of Projections 

A geometric approach will be taken to formulate the 
time-update for the least squares lattice algorithm. A 
simple example of projecting a vector |y> T onto a vector 
|x><p, in a three dimensional real space (Figure A-l) , by 
decomposing |x> T into its past and present components will 
facilitate the development of the projection time-update. 
Using the coordinate projector defined in Equation (A. 26), 
two new vectors, containing the past and present components 
of |x> T , will be defined as 

| x n >T s PtI x >T = [0,0, x t1 ' (A. 47) 

and 


| X - > 'p - | X>T — [XQ,Xp, 


• - • t 


x-p - 1 , o ] ' . 


(A. 48) 
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Several important relationships pertaining to Equations 
(A. 47) and (A. 48) are ^Yn^n^T = <Ynl x> T = < yl x n > T anc ^ 
<y. |x_>ip = <y|x_>T = <y.|x>-p, respectively. 

Taking a step in the opposite direction, the vector ly>T 
may be decomposed into ly n >T and IY->t and then projected 
onto |x>t such that 

P x |y> T = P x I y- >T + p xlYn > T • (A. 49) 

At this point the time index, T, will be dropped. Expanding 
Equation (A. 49) yields 

I x> <x | x> " l<x | y> = | x> <x | x> " l<x | y. > + | x> <x | x> * ^ <x | y n > 

= | x> <x | x> ' !<x. | y> + | x> <x | x> * 1 <x n | y> . 

(A. 50) 

The second equality of Equation (A. 50) reveals a 
decomposition of |x> into |x n > and |x_>. 

In order to develop a projection update in terms of | x n > 
and |x.>, an oblique projection needs to be defined. An 
oblique projection is a vector composed of the product of the 
magnitude of an orthogonal projection of |y> onto |x.> and 
I x> . The oblique projection takes the form 

P x - ly> B I x> <x. | X- > ■ 1<X- | y> . (A. 51) 

A correction component based on the present time sample |x n > 
is needed to complete the projection of |y> onto | x> . 

The decomposition of P x |y> based on the oblique 
projection in Equation (A. 51) yields 



p x iy> = p x -Iy> + p x PTP- L x-ly> 


(A. 52) 


where the second term on the right-hand side of Equation 
(A. 52) is the correction component based on the present time 
sample. A derivation of Equation (A. 52) is given as a 
geometric formulation in the next section. 

Geometric Formulation 

The first term on the right-hand side of Equation (A. 52) 
can be expressed as 

P x -ly> = IxXx. |x->" 1 <x. I y. > . (A. 53) 

Using previously defined relationships, the vector |y> in 
Equation (A. 51) has been replaced by the vector ly.> to form 
Equation (A. 53). The inner products in Equation (A. 53) 
represent the coefficient defined at time T-l. This 
relationship can be seen in Figure A-l where a is P x _ |y>. 

As stated previously the second term on the right-hand 
side of Equation (A. 52) is a correction factor. Referring to 
Figure A-l, the term P-*- x . |y> is the vector b where 

p_L x-ly> = ly> - p x - ly> • (a. 54) 

The coordinate projector P>p allows for the extraction of the 
most recent component. Therefore, P-L x _|y> is projected onto 
Py to form c. in Figure A-l. Finally, the result can be 
projected onto |x> to form P x Pij<P-L x . | y> . This correction 
factor d can now be added to a to form P x |y>. 
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At this point a more rigorous proof that Equation (A. 52) 
is equal to Equation (A. 50) is given. Taking the first term 
on the right-hand side of Equation (A. 50) 

|x><x|x> ' 1 <x|y> = |x><x|x>' 1 [<x. |x><x. lx)' 1 ] <x. |y> 

= | x> <x | x> " 1 [<x |x> - <x n | x> ] 

<x . | x> ” 1 <x - I y> 

= IxXx. |x->' 1 <x. |y> - | xXx | x> * 1 

<x n | xXx. |x_>‘ 1 <x. | y> 

= P x -ly> ‘ p x p T p x- ly> (A. 55) 

and adding the second term of Equation (A. 50) to (A. 55) gives 

p x iy> = p x -iy> • p x p T p x- iy> + p x p T I y> 

= p x -ly> + P x P T P J - x -|y> . (A. 56) 

From Equation (A. 56), the projection orthogonal to P x ly> 
can be written as 

p - L x iy> = pJ -x-ly> - p x p T pJ -x-ly> - (a. 57 ) 

The inner product of |y> with Equation (A. 57) will lead to 
the desired form needed for the reflection coefficient and 
error covariance time-updates. Therefore, taking the inner 
product and replacing the time indexes results in 

< ylE J - x ly>T = <yl p - L xly > T-l + <ylP J - x ln> T <n| P-*- x . |y> T . 

(A. 58) 

This can be proven by writing Equation (A. 57) as 
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p- L x iy> = p- L x-iy> - PT pl x-iy> + PT pJ -x-iy> - p x p TP- L x-iy> 

(A. 59) 

and using the definition of the orthogonal projection to 
yield 


p-L x ly> = p^tp^x- ly> p^"xPtP^"x- I y ^ (a.60) 

The result of the inner product |y> and Equation (A.60) can 
be written as 

<y|p j - x ly> = <yl P- L T pi 'x-ly > + <yl p j -x p t p ‘ L x- ly> • (a.6d 

The first term on the right hand side of Equation (A. 61) is 
the inner product of the error in projecting ly.> onto |x.> 
which is illustrated in Figure A-l. Given this projection, 
the time -update for the inner product can be written as 

<y i pJ -xly > T = <y l pJ -xly>T- 1 + <yl pJ_ xl n> T <n l p_L x- ly>T • 

(A. 62) 

Angle Between Two Subspaces 

The concept of an angle between two subspaces will be 
developed in this section in order to expand on the time- 
update of the inner product in Equation (A. 58) . Referring to 
Figure A-l let the angle between |x_> and |x> be denoted by 
0. The geometric relation can be expressed as 

cos^S = 1 - <x n |x n >/<x | x> = <x. | X. >/<x | x> (A. 63) 


or 
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sin 2 0 = <x n |x n >/<x|x> . (A. 64) 

An alternative form of sin 2 0 can be written in terms of 
projections defined in the previous section where 


sin 2 0 = <n|*P x |n> . (A. 65) 

Considering an n-dimensional complex space, Xi, n , Equation 
(A. 64) can be written as 


sin20 1(n = <n|X lfn ><X lfn |Xi >n >-l<X 1/n |n> 

= <n|*Pp /n |n> (A. 66) 

or Equation (A. 63) as 

cos 2 0p /n = <n | *P-^i t n I n> . (A. 67) 

An order-update for cos 2 0p )n ,T can be found using Equation 
(A. 37) where 


cos 2 0i f n +l , T = cos20 l,n,T ' r *n , T - 1 R ' r n , T - l r n , T - 1 . 

(A. 68) 

Another update that will be needed in developing an algorithm 
to determine the filter coefficients is 


cos 2 0q ( n,T = cos2 0O,n-l,T ' r *n, T R ' r n, T r n, T (A. 69) 

where it should be noted that cos 2 0p = cos 2 0g , n - 1 , T - 1 • 

Also, for notational purposes in later sections let y n< T = 
cos 2 0q , n- 1 , T • 

Having defined cos 2 0, it is possible to show the 
relationship between the orthogonal and oblique projections 
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on |x>. Multiplying Equation (A. 57) by P<p and using Equation 
(A. 67), the relationship is found to be 

p T pl xly> = p T pJ -x- ly> - PtPxPt^x- l x> 

= [Prp-P T P X ] P-pP-L- x . |y> 

= PtP-'-xPtP^x- ly> 

= PtP^x- I y>cos 2 0 . (A. 70) 

Exact Time -Update Formula 

Using Equations (A. 58) and (A. 70), the time update of 
the inner product using the orthogonal projection, P-Lj n <p 
can now be written as 

O’Cu | p ^"l , n , T I s ^u| P“^l , n, T I V ^T ” 

<u| *P^"i , n , T - 1 1 V> T - 1 
= <u| *P- L l, n , T ln> T 

<n | * p_L i , n , t I v>sec2 6l , n, T (A. 71) 

where |u>, |v> € Hrp. 

Time -Updates 

Using Equation (A. 71), the time-update for the 
reflection coefficients can be written as 

k n+l,T = k n+l , T - 1 + <x l pJ 'l,n,Tl n> T 

<n I p_L l,n,Tl z' n ' 1 x> T sec 2 0i /ri/T 
= k n+l , T - 1 + e *n,T r n,T- lS ec 2 0i , n , T • (A. 72) 

The conjugate placed on the forward error in Equation (A. 72) 
is due to the complex inner product defined in Equation 
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(A. 19) . This point is made due to the difference seen when 
comparing Equation (A. 72) to that in Lee [22] . The 
covariance terms for the backward and forward errors can also 
be updated to form 

Re n , T = Re n,T-l + e *n, T e n, T sec20 l , n, T (A. 73) 

and 

Rr n,T = Rr n, T - 1 + r *n, T r n, T sec20 O , n- 1 , T • (A. 74) 

Sliding Exponential Windows 

A window to reduce the effect of past data samples can 
be incorporated into the algorithm through the definition of 
the complex inner product. Referring to the complex inner 
product defined in Equation (A. 19), a new inner product can 
be now be defined as 

T 

<x|y>T = Z x*ta T ' fc yt 0<a<l . (A. 75) 

t=0 

This definition of the complex inner product can be used to 
transform the time-update of the reflection coefficient in 
Equation (A. 72) to 

k n+l,T = ^n+l , T - 1 + e *n,T r n,T-l sec20 l,n,T • (A. 76) 

The normalizations in the next section will eliminate the 
window a from every order of the lattice structure except the 


zeroth . 
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The Complex Unnormalized Recursive Least 
Squares Lattice Estimation Algorithm 

Table A- 1 contains the necessary equations and 
initializations needed to implement the complex unnormalized 
recursive least squares lattice estimation algorithm. For 
notational purposes in Table A- I, let Yn,T = cos ^®0,n-l,T • 
The resulting lattice structure is given in Figure A- 2. 

Normalized Ladder Recursions 

In this section, a normalization of the recursive least 
squares lattice algorithm will be used to reduce the number 
of equations needed in for the time and order updates in 
Table A- I. The normalizations in this section will lead to a 
form of the algorithm which consists of a time-update for the 
reflection coefficient and an order -update for the backward 
and forward prediction errors. 

Variance Normalization 

In this section, the forward and backward prediction 
errors and the reflection coefficients will be normalized by 
the square root of their respective covariances. The 
normalized prediction errors are defined as 

|v n >T s I e n^T^ e n I e n" > T _ (A. 77) 

liln>T * I r n- > T^ r n I r n- > T" ■ (A. 78) 

The reflection coefficient is normalized in a similar manner 


and takes the form 
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TABLE A- I 

Complex Least Squares Adaptive Lattice Algorithm 


Input parameters: 

N = maximum order of lattice 
y T = data sample at time T 
a = exponential weighting factor 
Do for T = 0 to T max 
e 0 , T = r 0,T = X T 

rS 0,T = Rr 0,T = aRe 0,T-l + YTY*T 
Y0,T = 1 

Do for n = 0 to {min(N,T)-l} 

k n+l,T = ^n+l,T-l + e *n, T r n, T- l/Yn, T- 1 

Yn+1 , T =Yn,T * r *n, T R ' r n, T r n, T 

kr n+l,T = k n+l , T R ' r n, T- 1 

e n+l , T = e n, T ' k * r n+l , T r n, T- 1 

Re n+1 , T = Re n, T * kr n+l , T k *n+1 , T 

ke n+l,T = R " e n,T k n+l,T 

r n+l , T = r n, T - 1 - ke n+l,T e n,T 

Rr n+1,T = Rr n, T - 1 * k *n+l , T ke n+1 , T 

Note: Division by zero where y = y, R r , R e / set 1/y = 0. 

Initialize the variables k, r, R e , R r , and y to zero. 
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r n+l,T * <v n |z' 1 n n > T 

= <e n |e n > T * 1 / 2 k n+lfT <z' 1 r n |z' 1 r n > T ' 1 / 2 • (A. 79) 

Having defined the normalizations in Equation (A. 77), (A. 78) , 

and (A. 79), a normalization of the order -updates for the 
backward and forward prediction errors can now be developed. 
A useful relationship involving the order -update of the 
forward prediction error covariance in Equation (A. 41) will 
be needed to express the order-update in terms of the 
normalized variables defined in Equations (A. 77), (A. 78) , and 

(A. 79). Multiplying both sides of Equation (A. 41) by 
<e n l e n > T yields 

<e n+l l e n+l>T <e nl e n>T = < e n I e n>T <e n I e n>T ’ < e n I 2 ' lr n>T 

< e n I e n>T< z " lr n I z ~ lr n>T" 1 
<z' 1 r n |e n > T {A. 80) 

or 

<e n+ll e n+l>T< e nl e n>T = <e n I e n >T< e n I e n>T U ' < e nl e n>T' 1/2 

k n+l,T< z ~ lr nl z ~ 1 rn>T' 1/2 

<e n I e n > T 1 ‘^ 2k *n+l , T 

<z' ir n | z‘ 1 r n > T ' 1 / 2 ] . (A. 81) 

Using the normalization of the reflection coefficients in 
Equation (A. 79), Equation (A. 81) can be written as 

<e n+l I e n+l>T <e n I e n > T = < e n I e n>T <e n I e n>T 

[i - r* n+ 1 , T r n+1/T ] . (a. 82 ) 

Dividing both sides of Equation (A. 82) by <e n | e n > T < e n I e n > T 
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and taking the square root yields 

<e n+i I e n+l>T 1/,2<e nl e n > T' 1/ ' 2 = H ' r *n+l , T r n+1 , 1 ^ 2 • 

(A. 83) 

Taking the order -update of the forward prediction error in 
Equation (A. 40) and dividing both sides by <e n +i I e n+l > t^' 2 
yields 


l e n+l>T< e n+l l e n+l>T' 1/2 = < e n+l I e n+l>T' 1/2 <e n I e n>T 1/2 

[|e n > T <e n |e n > T ' 1 / 2 - 

|z- 1 r n >T< z ' lr nl z ' lr n>T ‘ 1 
k *n+l,T <e nl e n > T' 1 ^ 2 ^ • 

(A. 84) 

Using the normalizations defined previously and Equation 
(A. 83), Equation (A. 84) takes the form 

l¥n+l>T = Ilv n >T - I z ' 1 Hn > T r *n+l,T] 

[1- r* n+ i >T r n+ i /T ] -1/2 (A. 85) 

The same type of normalization applies to the order -update 
for the backward error which can be written as 

IHn+l > T = Nz' 1 n n >T ' I ^n > T r n+l , t] [1* r *n+l , T r n +1 , T^ " 

(A. 86) 

Angle Normalization 

The next normalization is defined as 


l v n>T = I Yn>T sec ®l , n, T 


(A. 87) 
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and 


I z " lf ln > T * I z ' 1 Hn > T sec0 l,n,T • (A. 88) 

This normalization defines the forward and backward 
prediction errors in terms that are no longer observable, but 
the definition of the reflection coefficients will remain 
unchanged. 

At this point, a useful relationship, needed in the 
normalization of the time-update of the reflection 
coefficients, will be developed. Applying the inner product 
defined in Equation (A. 75) , the time-update for the forward 
prediction error covariance in Equation (A. 73) takes the form 

a <e nle n >T-l = < e nl e n>T * <e n 1 n> T <n I e n > T sec 2 0i ; n? T . 

(A. 89) 

Factoring <e n |e n >T from the right hand side of Equation 
(A. 89) yields 

a<e n l e n>T-l = <e n I e n>T t 1 * < e nl e n > T ’ 1 

<e n |n>-p<nIe n >i.sec 2 0i <n> T] • (A. 90) 

Dividing both sides of Equation (A. 90) by <e n |e n ><p and 
grouping terms gives 

a<e n |e n > T -i<e n |e n > T ' 1 = 1 - <e n |e n > T ' 1 / 2 e* n/T sec0i / n /T 

<e n l e n>T' 1 / 2 e n,T sec 0 l,n,T - 


(A. 91) 
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Using the normalization defined in Equation (A. 87), Equation 
(A. 91) takes the form 

a<e n I ©n>T- l< e n I e n>T" 1 = [1 ' v *n,T v n,Tl • (A. 92) 

Finally, dividing by a and taking the square root of Equation 
(A. 92) yields 

< e nle n >T-l 1/2 < e nlen>T' 1/2 = «' 1/2 [1 - v* n< T v Rf T ] i/ 2 . 

(A. 93) 

Using the time-update for the backward prediction error 
covariance in Equation (A. 74) , a similar development yields 

<z‘ 1 r n | z' 1 r n > T . i 1 / 2 <z' 1 r n | z'^-rni-r' 1 / 2 

= a"l/ 2 [l - <z'ln n |n>T<n|z~ln n >T’] • (A. 94) 

Using the time-update for the reflection coefficient in 
Equation (A. 76) and the normalization in Equation (A. 79) , the 
normalized time-update of the reflection coefficient can be 
written as 

r n+l,T = <e nl e n > T’ 1//2<e nl e n > T-l 1 ^ 2<e nl e n > T-l" 1 ' /2 
o^n+l , T - l<z - ir n I z ' 1 r n >T - 1 ' 1/2 
<z' 1 r n |z'lr n > T .i 1 / 2 <z' 1 r n |z' 1 r n > T ' 1 / 2 + 
<e n l e n>T' 1/2 < e nl n >T s ec0i,n,T 

<z" ^r n | z ' -*-r n >T" -*-/ 2 <n | z' lr n >>psec0i , n, T • 


(A. 95) 
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Equation (A. 95) can now be written in terms of normalized 
variables as 

r n+l,T = ' v *n,T v n,T] 1/,2r n+l,T-l 

t 1 ‘ n*n, T- lH n,T-l] 1,/2 + v *n,Tnn,T-l (A. 96) 

using Equations (A. 77), (A. 78), (A. 79), (A. 87), (A. 93), and 

(A. 94) . 

Ladder Recursions 

At this point, a second normalization of the order- 
updates for the forward and backward prediction errors is 
possible. This normalization is based on the terms defined 
in Equations (A. 87) and (A. 88). Equation (A. 85), the 
normalized order-update for the forward prediction, can be 
normalized in the manner given in Equation (A. 87). Equation 
(A. 66) can be written as 

|Vn+l>Tsec6 1>n+1(T /sec0 1(n+1>T = secQ li n , T /sec0 1/ n , T 

[|v n >T ' l z ' lll n > T 
r*n+l , t! t 1 ' r*n+i,T 
r n +l,T ]* 1/2 • (A. 97) 

Using the normalizations in Equations (A. 87) and (A. 88), 
Equation (A. 97) takes the form 

l v n+l>T sec0 l,n,T/ sec0 l,n+l,T = 1 1 v n>T * I z * lT1 n>T r *n+l , 

ti-rVi, T r n+ i,T ]- 1/2 • 


(A. 98) 
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In order to complete the normalization, a relationship 
between sec0^ t n -p/secO^ , n+1 , T an< 3 the normalized variables 
must be developed. Writing the order-update for cos 2 0i /rl/T 
in Equation (A. 68) as 

cos 2 0i, n +l,T/ cos2 01,n,T = t 1 ' r *n , T - l cos ' l0 l , n , T 

R r n,T- l r n,T- l co s’^e 1 , n , T ] 

(A. 99) 

will lead to the desired relationship. The trigonometric 
relationship 

sec 2 0i >ri( T = l/cos 2 0p ' n/ -j. (A. 100) 

is stated here as a reminder. Using the relationship stated 
in Equation (A. 100) and the normalizations given in Equations 
(A. 77), (A. 78) , and (A. 87), Equation (A. 99) takes the form 

sec 2 0i, n ,T/ sec20 l,n+l,T = t 1 ' H*n, T- l^n, T- ll • (A. 101) 

Using the relationship given in Equation (A. 101), the 
normalized order -update of the forward prediction error in 
Equation (A. 98) can be written as 

v n+l , T = t 1 ' r *n+l , T r n+1 , t! t v n,T ' r * n+1 , T^n, T - 1^ 

[1 - T^n.T-l^T-l] " 1 ^ 2 (A. 102) 

while looking only at the Tth component. A similar 
development leads to the normalized order -update of the 
backward prediction error which can be written as 
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Hn+1 , T = Cl ’ r *n+l,T r n+l,Tl ' 1/2 [t1n,T-l ' r n+l,T v n,T] 

[1 - V* n , T v n , t] " • (A. 103) 

The Complex Normalized Recursive Least 
Squares Lattice Estimation Algorithm 

The order -updates for the backward and forward 
prediction error covariances and for the trigonometric 
relationship cos^S have been embedded in the normalizations 
defined in this section. A reduction in the number of 
equations, from six to three, required in the lattice 
recursions has been achieved. Table A- II gives the complex 
square root normalized recursive least squares lattice 
estimation including the necessary initializations, and 
Figure A- 3 gives the resulting lattice structure. 

Complex Prediction Error Filter Coefficients 
Introduction 

An algorithm for determining the coefficients of the 
prediction error filter is given by Honig [18] as a lattice 
structure implementation for the real data case. The complex 
form of this algorithm is presented in this section to 
support the complex square root normalized recursive least 
square lattice estimation algorithm. A transfer function 
relationship will be defined for the forward and backward 
prediction errors which can then be implemented as a lattice 
structure. The coefficients of the filter can then be 


obtained from the transfer function. 
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TABLE A -II 

Complex Square Root Normalized Recursive Least 
Squares Lattice Estimation Algorithm 


Initialize : 

R-l = a a is a small positive value 

for T = 0,T max 

Rrp = CCR/p _ 2 + Y *pyip 

V 0,T = ^0,1 = Yt( r T^ 

for n = 0 to [min{T,N} -1] 

r n+l,T = t 1 ' * * * * v *n,T v n,T3 1 ^ 2 r n+l,T-l 

t 1 ' T l*n,T-l T ln,T-l] 1 ^ 2 + v *n,Tn n ,T-l 

v n+l , T = t 1 ' r* n +l , T^n+1 , t v n, T ' r *n+l , T%, T- 1^ 

t 1 ' n*n,T-l r ln / T-ll ~ 1 ^ 2 

’In+l, T = t 1 ' r*n+l t Tr n +i , ' '*'^ 2 t^n, T- 1 ' r n+l,T v n,T^ 

t 1 * v *n,T v n,T^ ‘ 1 ^ 2 

Note: Division by zero where y = 1/x : x = 0 should result in 
y = 0. Initialize the variables T, v, and n to zero. 
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Transfer Functions 

The nth order forward prediction error, at time T, can 
be written as 

e n, T = x T + a n , l , T x t . i + ... + a n , n , T x T . n . (A. 104) 

Taking the Z- transform of this forward prediction error and 

dividing through by x(z) , the transfer function can be 
defined as 

A n , T ( z ) s 1 + a n, 1 , T 2 ' 1 + ••• + a n,n,T z ~ n • (A. 105) 

A n,T^ z ^ can be written in vector form as 

A n , T ( z > = 1 + A'n,T z ' lz n (A. 106) 

where 


A n,T - t a n, 1 , T ••• a n,n,T 

and 

z n = [ 1 z'l ... z‘ n +l] ' . 

Similarly, the nth order backward prediction error, at time 
T, can be written as 

r n, T = x T-n + b n,0,T x T + + b n,n-l,T x T-n+l (A. 107) 

where the transfer function for this filter is 

B n,T(z) = z ' n + b n, 0 , T + ••• + b n, n- 1 , T z ' n+1 • (A. 108) 

The vector form of B R/T (z) can be written as 
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B n , T (z) = z' n + B 'n,T z n (A. 109) 

where 

B n , T = [ b n , 0 , T • * • b n , n - 1 , T 1 

A version of the backward prediction error based on 
estimating XT- n using the coefficients obtained at time T-l 
will be needed in the order updates of the transfer functions 
and can be written as 

r n, t| T- 1 = x T-n + b n,0,T-l x T 

+ ... + b n,n-l,T-l x T-n-l (A. 110) 

and the associated transfer function as 

B n ,T-l<z) = [2' n + B "n,T- lZ.n) (A. Ill) 

where 


B n,T-l ~ [ b n , 0 , T - 1 ••• b n,n-l,T-l 1' • 

An order update for the forward and backward prediction 
errors is given in Equations (A. 40) and (A. 44), respectively. 
The transfer functions for the forward error order update, 
e n+l,T' an< 3 tbe backward error order update, r n+ i t, can be 
written as 

A n+ 1 , T ( 2 ) = A n,T^ z ) ' z ’ lB n , T - 1 ( z ) R ’ r n , T - 1 k *n + 1 , T 

(A. 112) 


and 
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B n+1 , T ~ z ’ lB n,T-l (z ) ' A n, T (z) R " e n, T k n+1 , T (A. 113) 
respectively. 

In order to implement the transfer functions in terms of 
a lattice structure, a relationship between B n> -p(z) and 
B n,T-l( z ) needs to be developed. The time -update formula 
given in Equation (A. 57) can be used to form the time-update 

pJ "0,n-l,T = P-^O, n-l,T " PQ^-l^PTP-^O.n-l.T (A. 114) 

where P.-Lq n -l,T is the orthogonal oblique projection 
operator. Using Equation (A. 70), Equation (A. 114) can be 
written as 

P-^n-^T = P-^n-^T ‘ 

p 0 , n - 1 , T p T p_ ^0 , n - 1 , T sec ^®0 , n - 1 , T • (A. 115) 

The backward prediction error at time T, r n+l,T' can ke 
obtained by pre and post multiplying Equation (A. 115) by t < 7T | 
and | z' n x>T, respectively, to form 

T <7r | P_L 0 , n - 1 , T I z " nx> T = T <7r | p_ *"0 , n- 1 , T- 1| z' n x>T ' 

T<^|P0,n-l,T p TP 1 0,n-l,T |z' n x> T 
sec^Sg , n . , t • (A. 116) 

From Equation (A. 30), the term on the left hand side of the 
equals sign in Equation (A. 116) is the backward prediction 
error at time T. The first term on the left hand side of the 
equals sign represents the estimate of the backward 
prediction error at time T using the coefficients defined at 
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time T-l. The second term on the right side of Equation 
(A. 116) is not so apparent. The second term can be written 
as 

T <7r | POfn-lfTPT^Ojn-lfT I z' n x> T sec 2 0 O , n -l,T 

= T <7r | p 0,n- 1, t| 7r> T <7r | *T P ^'0 , n- 1 , t| z " nx> T 

sec 2 0o, n -l,T. (A. 117) 

From Equation (A. 30), Equation (A. 117) can be written as 

T <7r | p 0,n-l,T p T p_L 0,n-l,T I z ' n x> T sec 2 0Q , n - 1 , T 

= T <7r | p O,n-l,T| 7r> T r n,T sec20 O,n-l,T. (A. 118) 

The projection of | 7T> «j, on to the subspace Xq^.^t in 
Equation (A. 118) can be written as 

T <7r | p 0,n-l,T| ff> T = T <7r | | x 0,n-l>T <x 0,n-l| X 0,n-1 > T' 1 

<x 0 , n- 1 1 T* | ^T . (A. 119) 

Using the least squares estimate defined in Equation (A. 11), 
Equation (A. 119) can be written as 

T<H p 0 , n- 1 , t| ^>T = T <7r | | x 0,n-l > T c n,T (A. 120) 

where 

Cn,T = f c n,0,T ••• C n,n-1,T 1 
Equation (A. 118) can now be written as 
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T <7r | p O,n-l,T| 7r> T r n,T sec20 O,n-l,T = T <7r | I x 0 , n- l>T c n, T 

r n, T sec2 0 O , n- 1 , T . 

(A. 121) 

Taking the z- transform of Equation (A. 121) and dividing 
through by x(z) yields 

C n , T (z) = c n , 0 + c n, l z ' 1 + ... + c n, n- l z " n+ ^ 

= C n t T— n - 1 • (A. 122) 

The z-transform of Equation (A. 116) divided by x(z) can now 
be written as 

®n,T^ z ^ = ®n,T-l^z) - C n/ <p (z) r n , T/Yn, T (A. 123) 

where 


Y n ,T = cos 2 0 O/ n-l,T • 

An order -update of C n rp(z) is needed to implement the 
lattice structure. Using the concept of the decomposition of 
subspaces defined in Equation (A. 33), the projection operator 
Po,n,T can ke written as 

p 0,n,T = P 0,n-1,T + I r n > T <r n| r n > T’ 1<r n| T* • (A. 124) 

Pre and post multiplying Equation (A. 124) by t< 7 t| and | tt> T , 
respectively, yields 

T <7r | p 0 , n, t| 7T>T = T <7r | p 0 , n- 1, t| ^T + 

T <7r | | r n > T <r n| r n>T’ 1<r n| T*| ^T • (A. 125) 


Using the transfer functions defined in Equations (A. 108) and 
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(A. 122), the transfer function associated with Equation 
(A. 125) can be written as 

C n +l,T(z) = C n , T (z) + B n , T (z)R- r n/T r* n(T . (A. 126) 

Unnormalized Prediction Error 
Filter Coefficients 

The coefficients of the prediction error filter can be 
obtained from the lattice implementation of the transfer 
functions by associating terms in the transfer functions that 
are multiplied by the same power of the variable z‘P p = 
0,...,n. Table A-III gives the algorithm for the lattice 
implementation of the coefficients associated with the 
prediction error filter, and Figure A-4 gives the associated 
lattice structure. 

Normalized Transfer Functions 

A normalized lattice structure implementation of the 

transfer functions which involves the normalized forward and 
backward prediction errors and the normalized reflection 
coefficient defined for normalized lattice structure given in 
the previous section will be developed. The normalized 
transfer functions will be defined as 


A n , t ( 2 ) 

= A n>T (z) <e n | 

en>T - 1/2 

(A. 127) 

Bn, T (z) 

= B n / T ( z ) <r n | 

r„> T -l/2 

(A. 128) 


Bn,T-l< z > = Bn,T-l <r n| r n>T-l‘ 1/2 (A. 129) 


and 
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TABLE A- III 

Algorithm for Determining the AR Model Coefficeints 
Based on a Complex LS Adaptive Lattice Structure 


For i = 0,...,N do the following 
b p,-l,T = 0 for P = 0,...,N-1 
a 0 , 0 , T = b 0,0,T = 1 for i = 0 
c 0 , 0 , T = 0 for i = 0 

a 0 , i , T = b 0,i,T = c 0 , i , T = 0 for i > 0 
For p = 0 , . . . , N - 1 

b p,i,T-l = b p , i , T + c p, i,T r p,T/Yp,T 
c p+l,i,T = c p,i,T + b p,i,T R r p,T r p, T 
a p+l,i,T = a p,i,T ' b p,i-l,T-l R r p,T-l^ p+1 , T 
b p+l , i , T = b p, i - 1 , T- 1 ’ a p,i,T R ' e p,T k p+l,T 
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c n , T ( 2 ) “ c n , T ( 2 ) sec:0o , n - 1 , T . (A. 130) 

Based on Equations (A. 127), (A. 129), (A. 79), and (A. 82), 

Equation (A. 112) can be normalized as 

A n+ 1 , T ( 2 ) = [ A n , T ( 2 ) ' 2 ’ le n, T- 1 ^ 2 ^ r* n +i, t! 

[i-r*n+i,Tr n +i,T]' 1/2 • (a. 131 ) 

A similar normalization holds for Equations (A. 113) and where 

Bn+i,T( z ) = ' A n, T < z ) r n+l , t! 

[l-r* n +l,Tr n +l,Tl ' 1/2 • (A. 132) 

Using the previously defined normalizations and Equation 
(A. 94) with the time index T replaced with T+l, Equation 
(A. 123) can be normalized as 

Bn,T-l( z ) = [Bn,!^ 2 ) + c n , T < 2 > Hn, tJ 
' T l*n, T^n, t] ' 1 ^ 2 - 

The final normalization requires the relationship 

sec 2 0 o ,n-l,T/ sec2e O,n,T = t 1 *n*n,Tn n ,T] 
to form the normalized transfer function 

c n+l , T ( z ) = t c n , T ( 2 ^ + Bn, T ( z ) T l*n, 

[l-n*n, T n n ,T] ' 1/2 • (A. 135) 

The coefficients of the prediction error filter can be 
obtained from the normalized lattice implementation of the 


(A. 133) 


(A. 134) 


transfer functions in a similar manner as that used for the 



124 


unnormalized case. Table A-IV gives the algorithm for 
obtaining the filter coefficients from the normalized lattice 
structure and Figure A-5 gives the associated lattice 


structure . 
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TABLE A -IV 

Algorithm for Determining the AR Model Coefficients Based 
on a Normalized Complex LS Adaptive Lattice Structure 


Rrp — (XR'P . + X 'J'X'P 

For i = 0,...,N do the following 

b p,-l,T = 0 for p = 0 , . . . , N- 1 
a 0,0,T = b 0,0,T = Rt' 1 / 2 for i = 0 
C 0,0,T = 0 for i = 0 

a 0 , i , T = b 0 , i , T = c 0 , i , T = 0 for i > 0 
For p = 0, . . . ,N-1 

b p,i,T-l = [ b p,i,T + c p,i,Tnp,T^ 
t 1 ' , l*p / T T lp / T] " 1/,2 

c p+ 1 , i , T = t c p,i,T + b p,i,Tn*p,T] [ 1 -n*p,T T lp,T^ ' 1/2 
a p+ 1 , i , T = t a p, i , T ‘ b p, i- l,T-l r *p+l,T^ 

[i*r* p +i / T r P +i, t] " 1//2 

b p+l,i,T = [ b p,i-l,T-l ■ a p, i,Trp+l,T] 
t 1 ' r *p+l,T r p+l,T] ' 1 ^ 2 

where 

a p, i = a p, i/ a p, 0 



(b) A single section. 


Figure A- 5. Transfer function implemented as a 1 
structure for the complex square root normalized 
squares adaptive prediction error filter case. 


at tice 
least 











Appendix B 

NASA Model Parameters 


SIMULATION PARAMETERS 

A/C Distance to touchdown (kM) 
Aircraft Velocity (kts) 
Glideslope Angle (deg) 

Roll Attitude (deg) 

Pitch Attitude (deg) 

Yaw Attitude (deg) 

Az Integration Range/2 (deg) 

Az Integration Increment (deg) 

El Integration Range/2 (deg) 

El Integration Increment (deg) 
Rng Integration Increment (m) 
Random Number Seed (0-1) 

MICROBURST & CLUTTER 

Along Track Offset from TD (km) 
Cross Track Offset from TD (km) 
Rain Standard Deviation (m/s) 
Clutter Standard Deviation (m/s) 
Clutter Calc. Flag 
Reflectivity Calc. Thres . (dBz) 


Minimum Reflectivity (dBz) 
Attenuation Code (0,1,2) 

RADAR PARAMETERS 

Initial Radar Range (km) 

Number of Range Cells 
Antenna Az - if no scan (deg) 
Azimuth Scan Range/2 (deg) 
Azimuth Scan Increment (deg) 
Antenna Elevation (deg) 
Transmitted Power (watts) 
Frequency (GHz) 

Pulse Width (microsecs) 

Pulse Interval (microsecs) 
Receiver Noise Figure (dB) 
Receiver Losses (dB) 

Antenna Type (l=para., 2=flat) 


7.0 

150.0 

3.0 
0.0 
0.0 
0.0 

6.0 
0.3 
4.0 
0.2 

100.0 
0.224 


- 2.0 

0.0 

1.0 

0.5 

(l=on, O=off) 

1.0 (wet) 
-20.0 (dry) 
200.0 (clutter) 
-15.0 

2.0 


1.0 

40.0 

0.0 

0.0 

3.0 

1.0 

2000.0 

9.3 

1.0 

268.6 

4.0 

3.0 

1.0 
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Fourier Spectral Estimates of a Dry 
Microburst Plus Clutter 
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Fourier Spectral Estimates of a Dry 
Microburst Without Clutter 
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Fourier Spectral Estimates of Clutter 






















Appendix F 


Magnitude Response of the 10th Order 
FIR Filters 
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Phase Response of the 10th Order 
FIR Filters 
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Appendix H 


Filtered Spectrum Using the Appropriate 
Model Based Filter Coefficients 
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Filtered Spectrum Using the Filter 
Coefficients for Range Cell 20 
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Filtered Spectrum Using a Pulse 
Canceller Filter 
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